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 | |
---|
13 | %======================================================================= |
---|
14 | % Copyright 2008-2015, 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 | |
---|
31 | function DataOut=signal_FFTMean(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 = {'not used'};% 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 | |
---|
90 | %% get the variable to process |
---|
91 | Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse |
---|
92 | if isfield(Param.TransformInput,'IndexRange') |
---|
93 | IndexRange=Param.TransformInput.IndexRange; |
---|
94 | switch size(IndexRange,1) |
---|
95 | case 3 |
---|
96 | Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2),IndexRange(3,1):IndexRange(3,2)); |
---|
97 | case 2 |
---|
98 | Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2)); |
---|
99 | case 1 |
---|
100 | Var=Var(IndexRange(1,1):IndexRange(1,2)); |
---|
101 | end |
---|
102 | end |
---|
103 | np=size(Var);%dimensions of Var |
---|
104 | if ~isvector(Var) |
---|
105 | Var=reshape(Var,np(1),prod(np(2:end)));% reshape in a 2D matrix with time as first index |
---|
106 | end |
---|
107 | Var=Var-ones(np(1),1)*nanmean(Var,1); %substract mean value (excluding NaN) |
---|
108 | |
---|
109 | %% look for 'time' coordinate |
---|
110 | VarIndex=find(strcmp(Param.TransformInput.VariableName,DataIn.ListVarName)); |
---|
111 | TimeDimName=DataIn.VarDimName{VarIndex}{1}; |
---|
112 | TimeVarNameIndex=find(strcmp(TimeDimName,DataIn.ListVarName)); |
---|
113 | if isempty(TimeVarNameIndex) |
---|
114 | Time=1:np(1); |
---|
115 | TimeUnit='vector index'; |
---|
116 | else |
---|
117 | Time=DataIn.(DataIn.ListVarName{TimeVarNameIndex}); |
---|
118 | TimeUnit=['Unit of ' TimeDimName]; |
---|
119 | end |
---|
120 | % check time intervals |
---|
121 | diff_x=diff(Time); |
---|
122 | dx=min(diff_x); |
---|
123 | freq_max=1/(2*dx); |
---|
124 | check_interp=0; |
---|
125 | if diff_x>1.001*dx % non constant time interval |
---|
126 | check_interp=1; |
---|
127 | end |
---|
128 | |
---|
129 | %% claculate the spectrum |
---|
130 | specmean=0;% mean spectrum initialisation |
---|
131 | cospecmean=0; |
---|
132 | NbNan=0; |
---|
133 | NbPos=0; |
---|
134 | np_freq=floor(size(Var,1)/2); |
---|
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 |
---|
148 | sample(isnan(sample))=[]; |
---|
149 | end |
---|
150 | |
---|
151 | fourier=fft(sample);%take fft (complex) |
---|
152 | spec=abs(fourier).*abs(fourier);% take square of the modulus |
---|
153 | spec=spec(1:np_freq,:);%keep only the first half (the other is symmetric) |
---|
154 | specmean=spec+specmean; |
---|
155 | % end |
---|
156 | end |
---|
157 | specmean=specmean/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 | set(gca,'YLim',[1.0000e-06*max(specmean) 1.1*max(specmean)]) |
---|
169 | title (['power spectrum of ' Param.TransformInput.VariableName ]) |
---|
170 | xlabel(['frequency (cycles per ' TimeUnit ')']) |
---|
171 | ylabel('spectral intensity') |
---|
172 | legend({'spectrum','cospectrum t t-1'}) |
---|
173 | get(gca,'Unit') |
---|
174 | sum(specmean) |
---|
175 | if NbPos~=size(Var,2) |
---|
176 | disp([ 'warning: ' num2str(size(Var,2)-NbPos) ' NaN sampled removed']) |
---|
177 | end |
---|
178 | if NbNan~=0 |
---|
179 | disp([ 'warning: ' num2str(NbNan) ' NaN values replaced by linear interpolation']) |
---|
180 | %text(0.9, 0.5,[ 'warning: ' num2str(NbNan) ' NaN values removed']) |
---|
181 | end |
---|
182 | grid on |
---|
183 | |
---|
184 | |
---|