[753] | 1 | % 'signal_spectrum': calculate and display spectrum of the current field |
---|
| 2 | % operate on a 1D signal or the first dimension of a higher dimensional matrix (then average over other dimensions) |
---|
| 3 | % this function aplies the Welch method and call the function of the matlab signal processing toolbox |
---|
| 4 | % |
---|
| 5 | % OUTPUT: |
---|
| 6 | % DataOut: if DataIn.Action.RUN=0 (introducing parameters): Matlab structure containing the parameters |
---|
| 7 | % else transformed field, here not modified (the function just produces a plot on an independent fig) |
---|
| 8 | % |
---|
| 9 | % INPUT: |
---|
| 10 | % DataIn: Matlab structure containing the input field from the GUI uvmat, DataIn.Action.RUN=0 to set input parameters. |
---|
| 11 | % Param: structure containing processing parameters, created when DataIn.Action.RUN=0 at the first use of the transform fct |
---|
| 12 | |
---|
[810] | 13 | %======================================================================= |
---|
| 14 | % Copyright 2008-2014, LEGI UMR 5519 / CNRS UJF G-INP, Grenoble, France |
---|
| 15 | % http://www.legi.grenoble-inp.fr |
---|
| 16 | % Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr |
---|
| 17 | % |
---|
| 18 | % This file is part of the toolbox UVMAT. |
---|
| 19 | % |
---|
| 20 | % UVMAT is free software; you can redistribute it and/or modify |
---|
| 21 | % it under the terms of the GNU General Public License as published |
---|
| 22 | % by the Free Software Foundation; either version 2 of the license, |
---|
| 23 | % or (at your option) any later version. |
---|
| 24 | % |
---|
| 25 | % UVMAT is distributed in the hope that it will be useful, |
---|
| 26 | % but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
| 27 | % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
| 28 | % GNU General Public License (see LICENSE.txt) for more details. |
---|
| 29 | %======================================================================= |
---|
| 30 | |
---|
[753] | 31 | function DataOut=signal_spectrum(DataIn,Param) |
---|
| 32 | |
---|
| 33 | %% request input parameters |
---|
| 34 | if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0) |
---|
| 35 | VarNbDim=cellfun('length',DataIn.VarDimName); |
---|
| 36 | [tild,rank]=sort(VarNbDim,2,'descend');% sort the list of input variables, putting the ones with higher dimensionality first |
---|
| 37 | ListVarName=DataIn.ListVarName(rank); |
---|
| 38 | VarDimName=DataIn.VarDimName(rank); |
---|
| 39 | InitialValue=1;%default choice |
---|
| 40 | if isfield(Param,'TransformInput') && isfield(Param.TransformInput,'VariableName') |
---|
| 41 | val=find(strcmp(Param.TransformInput.VariableName,ListVarName)); |
---|
| 42 | if ~isempty(val); |
---|
| 43 | InitialValue=val; |
---|
| 44 | end |
---|
| 45 | end |
---|
| 46 | [s,OK] = listdlg('PromptString','Select the variable to process:',... |
---|
| 47 | 'SelectionMode','single','InitialValue',InitialValue,... |
---|
| 48 | 'ListString',ListVarName); |
---|
| 49 | if OK==1 |
---|
| 50 | VarName=ListVarName{s}; |
---|
| 51 | DataOut.TransformInput.VariableName=VarName; |
---|
| 52 | dlg_title = [mfilename ' calulates spectra along first dim ' VarDimName{s}{1}];% title of the input dialog fig |
---|
| 53 | prompt = {'nbre of points for the sliding window'};% titles of the edit boxes |
---|
| 54 | %default input: |
---|
| 55 | def={'512'};% window length |
---|
| 56 | np=size(DataIn.(VarName)); |
---|
| 57 | for idim=1:numel(np) % size restriction |
---|
| 58 | if idim==1 |
---|
| 59 | prompt=[prompt;{['index range for spectral dim ' VarDimName{s}{idim}]}];% titles of the edit boxes |
---|
| 60 | else |
---|
| 61 | prompt=[prompt;{['index range for ' VarDimName{s}{idim}]}];% titles of the edit boxes |
---|
| 62 | end |
---|
| 63 | def=[def;{num2str([1 np(idim)])}]; |
---|
| 64 | end |
---|
| 65 | if isfield(Param,'TransformInput') |
---|
| 66 | if isfield(Param.TransformInput,'WindowLength') |
---|
| 67 | def{1}=num2str(Param.TransformInput.WindowLength); |
---|
| 68 | end |
---|
| 69 | if isfield(Param.TransformInput,'IndexRange') |
---|
| 70 | for ilist=1:min(numel(np),size(Param.TransformInput.IndexRange,1)) |
---|
| 71 | def{ilist+1}=num2str(Param.TransformInput.IndexRange(ilist,:)); |
---|
| 72 | end |
---|
| 73 | end |
---|
| 74 | end |
---|
| 75 | num_lines= 1;%numel(prompt); |
---|
| 76 | % open the dialog fig |
---|
| 77 | answer = inputdlg(prompt,dlg_title,num_lines,def); |
---|
| 78 | DataOut.TransformInput.WindowLength=str2num(answer{1}); |
---|
| 79 | for ilist=1:numel(answer)-1 |
---|
| 80 | DataOut.TransformInput.IndexRange(ilist,1:2)=str2num(answer{ilist+1}); |
---|
| 81 | end |
---|
| 82 | end |
---|
| 83 | return |
---|
| 84 | end |
---|
| 85 | |
---|
| 86 | %% retrieve parameters |
---|
| 87 | DataOut=DataIn; |
---|
| 88 | WindowLength=Param.TransformInput.WindowLength; |
---|
| 89 | Shift=round(WindowLength/2);% shift between two windowsof analysis (half window length by default) |
---|
| 90 | |
---|
| 91 | %% get the variable to process |
---|
| 92 | Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse |
---|
[759] | 93 | if isfield(Param.TransformInput,'IndexRange') |
---|
| 94 | IndexRange=Param.TransformInput.IndexRange; |
---|
| 95 | switch size(IndexRange,1) |
---|
| 96 | case 3 |
---|
| 97 | Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2),IndexRange(3,1):IndexRange(3,2)); |
---|
| 98 | case 2 |
---|
| 99 | Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2)); |
---|
| 100 | case 1 |
---|
| 101 | Var=Var(IndexRange(1,1):IndexRange(1,2)); |
---|
| 102 | end |
---|
| 103 | end |
---|
[753] | 104 | np=size(Var);%dimensions of Var |
---|
| 105 | if ~isvector(Var) |
---|
| 106 | Var=reshape(Var,np(1),prod(np(2:end)));% reshape in a 2D matrix with time as first index |
---|
| 107 | end |
---|
[759] | 108 | Var=Var-ones(np(1),1)*nanmean(Var,1); %substract mean value (excluding NaN) |
---|
[753] | 109 | |
---|
| 110 | %% look for 'time' coordinate |
---|
| 111 | VarIndex=find(strcmp(Param.TransformInput.VariableName,DataIn.ListVarName)); |
---|
| 112 | TimeDimName=DataIn.VarDimName{VarIndex}{1}; |
---|
| 113 | TimeVarNameIndex=find(strcmp(TimeDimName,DataIn.ListVarName)); |
---|
| 114 | if isempty(TimeVarNameIndex) |
---|
| 115 | Time=1:np(1); |
---|
| 116 | TimeUnit='vector index'; |
---|
| 117 | else |
---|
| 118 | Time=DataIn.(DataIn.ListVarName{TimeVarNameIndex}); |
---|
| 119 | TimeUnit=['Unit of ' TimeDimName]; |
---|
| 120 | end |
---|
| 121 | % check time intervals |
---|
| 122 | diff_x=diff(Time); |
---|
| 123 | dx=min(diff_x); |
---|
| 124 | freq_max=1/(2*dx); |
---|
| 125 | check_interp=0; |
---|
| 126 | if diff_x>1.001*dx % non constant time interval |
---|
| 127 | check_interp=1; |
---|
| 128 | end |
---|
| 129 | |
---|
| 130 | %% claculate the spectrum |
---|
| 131 | specmean=0;% mean spectrum initialisation |
---|
| 132 | cospecmean=0; |
---|
| 133 | NbNan=0; |
---|
| 134 | NbPos=0; |
---|
| 135 | for pos=1:size(Var,2) |
---|
| 136 | sample=Var(:,pos);%extract sample to analyse |
---|
| 137 | ind_bad=find(isnan(sample)); |
---|
| 138 | ind_good=find(~isnan(sample)); |
---|
| 139 | if numel(ind_good)>WindowLength |
---|
| 140 | NbPos=NbPos+1; |
---|
| 141 | if ~isempty(ind_bad) |
---|
| 142 | sample=sample(ind_good); % keep only non NaN data |
---|
| 143 | NbNan=NbNan+numel(ind_bad); |
---|
| 144 | end |
---|
| 145 | %interpolate if needed |
---|
| 146 | if ~isempty(ind_bad)||check_interp |
---|
| 147 | sample=interp1(Time(ind_good),sample,(Time(1):dx:Time(end))); %interpolated func |
---|
[759] | 148 | sample(isnan(sample))=[]; |
---|
[753] | 149 | end |
---|
| 150 | spec=pwelch(sample,WindowLength);% calculate spectrum with Welch method |
---|
| 151 | cospec=cpsd(sample,circshift(sample,[1 0]),WindowLength);% calculate the cospectrum with the sample shifted by 1 time unit |
---|
| 152 | specmean=spec+specmean; |
---|
| 153 | cospecmean=cospec+cospecmean; |
---|
| 154 | end |
---|
| 155 | end |
---|
| 156 | specmean=specmean/NbPos; |
---|
| 157 | cospecmean=cospecmean/NbPos; |
---|
| 158 | |
---|
| 159 | %plot spectrum in log log |
---|
| 160 | hfig=findobj('Tag','fig_spectrum'); |
---|
| 161 | if isempty(hfig)% create spectruim figure if it does not exist |
---|
| 162 | hfig=figure; |
---|
| 163 | set(hfig,'Tag','fig_spectrum'); |
---|
| 164 | else |
---|
| 165 | figure(hfig) |
---|
| 166 | end |
---|
| 167 | loglog(freq_max*(1:length(specmean))/length(specmean),specmean) |
---|
| 168 | hold on |
---|
| 169 | loglog(freq_max*(1:length(cospecmean))/length(cospecmean),cospecmean,'r') |
---|
| 170 | hold off |
---|
| 171 | title (['power spectrum of ' Param.TransformInput.VariableName ]) |
---|
| 172 | xlabel(['frequency (cycles per ' TimeUnit ')']) |
---|
| 173 | ylabel('spectral intensity') |
---|
| 174 | legend({'spectrum','cospectrum t t-1'}) |
---|
| 175 | get(gca,'Unit') |
---|
[759] | 176 | sum(specmean) |
---|
| 177 | sum(cospecmean) |
---|
[753] | 178 | if NbPos~=size(Var,2) |
---|
| 179 | disp([ 'warning: ' num2str(size(Var,2)-NbPos) ' NaN sampled removed']) |
---|
| 180 | end |
---|
| 181 | if NbNan~=0 |
---|
| 182 | disp([ 'warning: ' num2str(NbNan) ' NaN values replaced by linear interpolation']) |
---|
| 183 | %text(0.9, 0.5,[ 'warning: ' num2str(NbNan) ' NaN values removed']) |
---|
| 184 | end |
---|
| 185 | grid on |
---|
| 186 | |
---|
| 187 | |
---|