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 | function DataOut=signal_FFTMean(DataIn,Param) |
---|
14 | |
---|
15 | %% request input parameters |
---|
16 | if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0) |
---|
17 | VarNbDim=cellfun('length',DataIn.VarDimName); |
---|
18 | [tild,rank]=sort(VarNbDim,2,'descend');% sort the list of input variables, putting the ones with higher dimensionality first |
---|
19 | ListVarName=DataIn.ListVarName(rank); |
---|
20 | VarDimName=DataIn.VarDimName(rank); |
---|
21 | InitialValue=1;%default choice |
---|
22 | if isfield(Param,'TransformInput') && isfield(Param.TransformInput,'VariableName') |
---|
23 | val=find(strcmp(Param.TransformInput.VariableName,ListVarName)); |
---|
24 | if ~isempty(val); |
---|
25 | InitialValue=val; |
---|
26 | end |
---|
27 | end |
---|
28 | [s,OK] = listdlg('PromptString','Select the variable to process:',... |
---|
29 | 'SelectionMode','single','InitialValue',InitialValue,... |
---|
30 | 'ListString',ListVarName); |
---|
31 | if OK==1 |
---|
32 | VarName=ListVarName{s}; |
---|
33 | DataOut.TransformInput.VariableName=VarName; |
---|
34 | dlg_title = [mfilename ' calulates spectra along first dim ' VarDimName{s}{1}];% title of the input dialog fig |
---|
35 | prompt = {'not used'};% titles of the edit boxes |
---|
36 | %default input: |
---|
37 | def={'512'};% window length |
---|
38 | np=size(DataIn.(VarName)); |
---|
39 | for idim=1:numel(np) % size restriction |
---|
40 | if idim==1 |
---|
41 | prompt=[prompt;{['index range for spectral dim ' VarDimName{s}{idim}]}];% titles of the edit boxes |
---|
42 | else |
---|
43 | prompt=[prompt;{['index range for ' VarDimName{s}{idim}]}];% titles of the edit boxes |
---|
44 | end |
---|
45 | def=[def;{num2str([1 np(idim)])}]; |
---|
46 | end |
---|
47 | if isfield(Param,'TransformInput') |
---|
48 | if isfield(Param.TransformInput,'WindowLength') |
---|
49 | def{1}=num2str(Param.TransformInput.WindowLength); |
---|
50 | end |
---|
51 | if isfield(Param.TransformInput,'IndexRange') |
---|
52 | for ilist=1:min(numel(np),size(Param.TransformInput.IndexRange,1)) |
---|
53 | def{ilist+1}=num2str(Param.TransformInput.IndexRange(ilist,:)); |
---|
54 | end |
---|
55 | end |
---|
56 | end |
---|
57 | num_lines= 1;%numel(prompt); |
---|
58 | % open the dialog fig |
---|
59 | answer = inputdlg(prompt,dlg_title,num_lines,def); |
---|
60 | DataOut.TransformInput.WindowLength=str2num(answer{1}); |
---|
61 | for ilist=1:numel(answer)-1 |
---|
62 | DataOut.TransformInput.IndexRange(ilist,1:2)=str2num(answer{ilist+1}); |
---|
63 | end |
---|
64 | end |
---|
65 | return |
---|
66 | end |
---|
67 | |
---|
68 | %% retrieve parameters |
---|
69 | DataOut=DataIn; |
---|
70 | WindowLength=Param.TransformInput.WindowLength; |
---|
71 | |
---|
72 | %% get the variable to process |
---|
73 | Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse |
---|
74 | if isfield(Param.TransformInput,'IndexRange') |
---|
75 | IndexRange=Param.TransformInput.IndexRange; |
---|
76 | switch size(IndexRange,1) |
---|
77 | case 3 |
---|
78 | Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2),IndexRange(3,1):IndexRange(3,2)); |
---|
79 | case 2 |
---|
80 | Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2)); |
---|
81 | case 1 |
---|
82 | Var=Var(IndexRange(1,1):IndexRange(1,2)); |
---|
83 | end |
---|
84 | end |
---|
85 | np=size(Var);%dimensions of Var |
---|
86 | if ~isvector(Var) |
---|
87 | Var=reshape(Var,np(1),prod(np(2:end)));% reshape in a 2D matrix with time as first index |
---|
88 | end |
---|
89 | Var=Var-ones(np(1),1)*nanmean(Var,1); %substract mean value (excluding NaN) |
---|
90 | |
---|
91 | %% look for 'time' coordinate |
---|
92 | VarIndex=find(strcmp(Param.TransformInput.VariableName,DataIn.ListVarName)); |
---|
93 | TimeDimName=DataIn.VarDimName{VarIndex}{1}; |
---|
94 | TimeVarNameIndex=find(strcmp(TimeDimName,DataIn.ListVarName)); |
---|
95 | if isempty(TimeVarNameIndex) |
---|
96 | Time=1:np(1); |
---|
97 | TimeUnit='vector index'; |
---|
98 | else |
---|
99 | Time=DataIn.(DataIn.ListVarName{TimeVarNameIndex}); |
---|
100 | TimeUnit=['Unit of ' TimeDimName]; |
---|
101 | end |
---|
102 | % check time intervals |
---|
103 | diff_x=diff(Time); |
---|
104 | dx=min(diff_x); |
---|
105 | freq_max=1/(2*dx); |
---|
106 | check_interp=0; |
---|
107 | if diff_x>1.001*dx % non constant time interval |
---|
108 | check_interp=1; |
---|
109 | end |
---|
110 | |
---|
111 | %% claculate the spectrum |
---|
112 | specmean=0;% mean spectrum initialisation |
---|
113 | cospecmean=0; |
---|
114 | NbNan=0; |
---|
115 | NbPos=0; |
---|
116 | np_freq=floor(size(Var,1)/2); |
---|
117 | for pos=1:size(Var,2) |
---|
118 | sample=Var(:,pos);%extract sample to analyse |
---|
119 | ind_bad=find(isnan(sample)); |
---|
120 | ind_good=find(~isnan(sample)); |
---|
121 | % if numel(ind_good)>WindowLength |
---|
122 | NbPos=NbPos+1; |
---|
123 | if ~isempty(ind_bad) |
---|
124 | sample=sample(ind_good); % keep only non NaN data |
---|
125 | NbNan=NbNan+numel(ind_bad); |
---|
126 | end |
---|
127 | %interpolate if needed |
---|
128 | if ~isempty(ind_bad)||check_interp |
---|
129 | sample=interp1(Time(ind_good),sample,(Time(1):dx:Time(end))); %interpolated func |
---|
130 | sample(isnan(sample))=[]; |
---|
131 | end |
---|
132 | |
---|
133 | fourier=fft(sample);%take fft (complex) |
---|
134 | spec=abs(fourier).*abs(fourier);% take square of the modulus |
---|
135 | spec=spec(1:np_freq,:);%keep only the first half (the other is symmetric) |
---|
136 | specmean=spec+specmean; |
---|
137 | % end |
---|
138 | end |
---|
139 | specmean=specmean/NbPos; |
---|
140 | |
---|
141 | %plot spectrum in log log |
---|
142 | hfig=findobj('Tag','fig_spectrum'); |
---|
143 | if isempty(hfig)% create spectruim figure if it does not exist |
---|
144 | hfig=figure; |
---|
145 | set(hfig,'Tag','fig_spectrum'); |
---|
146 | else |
---|
147 | figure(hfig) |
---|
148 | end |
---|
149 | loglog(freq_max*(1:length(specmean))/length(specmean),specmean) |
---|
150 | set(gca,'YLim',[1.0000e-06*max(specmean) 1.1*max(specmean)]) |
---|
151 | title (['power spectrum of ' Param.TransformInput.VariableName ]) |
---|
152 | xlabel(['frequency (cycles per ' TimeUnit ')']) |
---|
153 | ylabel('spectral intensity') |
---|
154 | legend({'spectrum','cospectrum t t-1'}) |
---|
155 | get(gca,'Unit') |
---|
156 | sum(specmean) |
---|
157 | if NbPos~=size(Var,2) |
---|
158 | disp([ 'warning: ' num2str(size(Var,2)-NbPos) ' NaN sampled removed']) |
---|
159 | end |
---|
160 | if NbNan~=0 |
---|
161 | disp([ 'warning: ' num2str(NbNan) ' NaN values replaced by linear interpolation']) |
---|
162 | %text(0.9, 0.5,[ 'warning: ' num2str(NbNan) ' NaN values removed']) |
---|
163 | end |
---|
164 | grid on |
---|
165 | |
---|
166 | |
---|