[876] | 1 | % 'signal_spectrum': calculate and display power spectrum of the current field |
---|
[753] | 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 | %======================================================================= |
---|
[924] | 14 | % Copyright 2008-2016, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France |
---|
[810] | 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 |
---|
[944] | 82 | if DataOut.TransformInput.IndexRange(1,2)-DataOut.TransformInput.IndexRange(1,1)<DataOut.TransformInput.WindowLength |
---|
| 83 | msgbox_uvmat('ERROR','WindowLength must be smaller than the total time index range') |
---|
| 84 | return |
---|
| 85 | end |
---|
| 86 | huvmat=findobj(allchild(0),'Tag','uvmat'); |
---|
| 87 | UvData=get(huvmat,'UserData'); |
---|
| 88 | Data=UvData.PlotAxes; |
---|
| 89 | YName=Data.ListVarName{1}; |
---|
| 90 | XName=Data.ListVarName{2}; |
---|
| 91 | yindex=DataOut.TransformInput.IndexRange(2,:); |
---|
| 92 | y=Data.(YName)(yindex); |
---|
| 93 | xindex=DataOut.TransformInput.IndexRange(3,:); |
---|
| 94 | x=Data.(XName)(xindex); |
---|
| 95 | haxes=findobj(huvmat,'Tag','PlotAxes'); |
---|
| 96 | axes(haxes); |
---|
| 97 | hbounds=findobj(haxes,'Tag','Bounds'); |
---|
| 98 | if isempty(hbounds) |
---|
| 99 | hbounds=rectangle('Position',[x(1) y(1) x(2)-x(1) y(2)-y(1)],'Tag','Bounds'); |
---|
| 100 | else |
---|
| 101 | set(hbounds,'Position',[x(1) y(1) x(2)-x(1) y(2)-y(1)]) |
---|
| 102 | end |
---|
[753] | 103 | end |
---|
| 104 | return |
---|
| 105 | end |
---|
| 106 | |
---|
| 107 | %% retrieve parameters |
---|
| 108 | DataOut=DataIn; |
---|
| 109 | WindowLength=Param.TransformInput.WindowLength; |
---|
| 110 | Shift=round(WindowLength/2);% shift between two windowsof analysis (half window length by default) |
---|
| 111 | |
---|
| 112 | %% get the variable to process |
---|
[944] | 113 | if ~isfield(DataIn,Param.TransformInput.VariableName) |
---|
| 114 | return |
---|
| 115 | end |
---|
[753] | 116 | Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse |
---|
[759] | 117 | if isfield(Param.TransformInput,'IndexRange') |
---|
| 118 | IndexRange=Param.TransformInput.IndexRange; |
---|
| 119 | switch size(IndexRange,1) |
---|
| 120 | case 3 |
---|
| 121 | Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2),IndexRange(3,1):IndexRange(3,2)); |
---|
| 122 | case 2 |
---|
| 123 | Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2)); |
---|
| 124 | case 1 |
---|
| 125 | Var=Var(IndexRange(1,1):IndexRange(1,2)); |
---|
| 126 | end |
---|
| 127 | end |
---|
[753] | 128 | np=size(Var);%dimensions of Var |
---|
| 129 | if ~isvector(Var) |
---|
| 130 | Var=reshape(Var,np(1),prod(np(2:end)));% reshape in a 2D matrix with time as first index |
---|
| 131 | end |
---|
[759] | 132 | Var=Var-ones(np(1),1)*nanmean(Var,1); %substract mean value (excluding NaN) |
---|
[753] | 133 | |
---|
| 134 | %% look for 'time' coordinate |
---|
| 135 | VarIndex=find(strcmp(Param.TransformInput.VariableName,DataIn.ListVarName)); |
---|
| 136 | TimeDimName=DataIn.VarDimName{VarIndex}{1}; |
---|
| 137 | TimeVarNameIndex=find(strcmp(TimeDimName,DataIn.ListVarName)); |
---|
| 138 | if isempty(TimeVarNameIndex) |
---|
| 139 | Time=1:np(1); |
---|
| 140 | TimeUnit='vector index'; |
---|
| 141 | else |
---|
| 142 | Time=DataIn.(DataIn.ListVarName{TimeVarNameIndex}); |
---|
| 143 | TimeUnit=['Unit of ' TimeDimName]; |
---|
| 144 | end |
---|
| 145 | % check time intervals |
---|
| 146 | diff_x=diff(Time); |
---|
| 147 | dx=min(diff_x); |
---|
| 148 | freq_max=1/(2*dx); |
---|
| 149 | check_interp=0; |
---|
| 150 | if diff_x>1.001*dx % non constant time interval |
---|
| 151 | check_interp=1; |
---|
| 152 | end |
---|
| 153 | |
---|
[944] | 154 | %% calculate the spectrum |
---|
[753] | 155 | specmean=0;% mean spectrum initialisation |
---|
| 156 | cospecmean=0; |
---|
| 157 | NbNan=0; |
---|
| 158 | NbPos=0; |
---|
| 159 | for pos=1:size(Var,2) |
---|
| 160 | sample=Var(:,pos);%extract sample to analyse |
---|
| 161 | ind_bad=find(isnan(sample)); |
---|
| 162 | ind_good=find(~isnan(sample)); |
---|
| 163 | if numel(ind_good)>WindowLength |
---|
| 164 | NbPos=NbPos+1; |
---|
| 165 | if ~isempty(ind_bad) |
---|
| 166 | sample=sample(ind_good); % keep only non NaN data |
---|
| 167 | NbNan=NbNan+numel(ind_bad); |
---|
| 168 | end |
---|
| 169 | %interpolate if needed |
---|
| 170 | if ~isempty(ind_bad)||check_interp |
---|
| 171 | sample=interp1(Time(ind_good),sample,(Time(1):dx:Time(end))); %interpolated func |
---|
[759] | 172 | sample(isnan(sample))=[]; |
---|
[753] | 173 | end |
---|
| 174 | spec=pwelch(sample,WindowLength);% calculate spectrum with Welch method |
---|
| 175 | cospec=cpsd(sample,circshift(sample,[1 0]),WindowLength);% calculate the cospectrum with the sample shifted by 1 time unit |
---|
| 176 | specmean=spec+specmean; |
---|
| 177 | cospecmean=cospec+cospecmean; |
---|
| 178 | end |
---|
| 179 | end |
---|
| 180 | specmean=specmean/NbPos; |
---|
| 181 | cospecmean=cospecmean/NbPos; |
---|
| 182 | |
---|
| 183 | %plot spectrum in log log |
---|
| 184 | hfig=findobj('Tag','fig_spectrum'); |
---|
| 185 | if isempty(hfig)% create spectruim figure if it does not exist |
---|
| 186 | hfig=figure; |
---|
| 187 | set(hfig,'Tag','fig_spectrum'); |
---|
| 188 | else |
---|
| 189 | figure(hfig) |
---|
| 190 | end |
---|
| 191 | loglog(freq_max*(1:length(specmean))/length(specmean),specmean) |
---|
| 192 | hold on |
---|
| 193 | loglog(freq_max*(1:length(cospecmean))/length(cospecmean),cospecmean,'r') |
---|
| 194 | hold off |
---|
| 195 | title (['power spectrum of ' Param.TransformInput.VariableName ]) |
---|
| 196 | xlabel(['frequency (cycles per ' TimeUnit ')']) |
---|
| 197 | ylabel('spectral intensity') |
---|
| 198 | legend({'spectrum','cospectrum t t-1'}) |
---|
| 199 | get(gca,'Unit') |
---|
[759] | 200 | sum(specmean) |
---|
| 201 | sum(cospecmean) |
---|
[753] | 202 | if NbPos~=size(Var,2) |
---|
| 203 | disp([ 'warning: ' num2str(size(Var,2)-NbPos) ' NaN sampled removed']) |
---|
| 204 | end |
---|
| 205 | if NbNan~=0 |
---|
| 206 | disp([ 'warning: ' num2str(NbNan) ' NaN values replaced by linear interpolation']) |
---|
| 207 | %text(0.9, 0.5,[ 'warning: ' num2str(NbNan) ' NaN values removed']) |
---|
| 208 | end |
---|
| 209 | grid on |
---|
| 210 | |
---|
| 211 | |
---|