source: trunk/src/transform_field/signal_FFTMean.m @ 766

Last change on this file since 766 was 766, checked in by sommeria, 10 years ago

sitgnal_FFTMean added

File size: 6.3 KB
Line 
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
13function DataOut=signal_FFTMean(DataIn,Param)
14
15%% request input parameters
16if 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
66end
67
68%% retrieve parameters
69DataOut=DataIn;
70WindowLength=Param.TransformInput.WindowLength;
71
72%% get the variable to process
73Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse
74if 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
84end
85np=size(Var);%dimensions of Var
86if ~isvector(Var)
87    Var=reshape(Var,np(1),prod(np(2:end)));% reshape in a 2D matrix with time as first index
88end
89Var=Var-ones(np(1),1)*nanmean(Var,1); %substract mean value (excluding NaN)
90
91%% look for 'time' coordinate
92VarIndex=find(strcmp(Param.TransformInput.VariableName,DataIn.ListVarName));
93TimeDimName=DataIn.VarDimName{VarIndex}{1};
94TimeVarNameIndex=find(strcmp(TimeDimName,DataIn.ListVarName));
95if isempty(TimeVarNameIndex)
96    Time=1:np(1);
97    TimeUnit='vector index';
98else
99    Time=DataIn.(DataIn.ListVarName{TimeVarNameIndex});
100    TimeUnit=['Unit of ' TimeDimName];
101end
102% check time intervals
103diff_x=diff(Time);
104dx=min(diff_x);
105freq_max=1/(2*dx);
106check_interp=0;
107if diff_x>1.001*dx % non constant time interval
108    check_interp=1;
109end
110
111%% claculate the spectrum
112specmean=0;% mean spectrum initialisation
113cospecmean=0;
114NbNan=0;
115NbPos=0;
116np_freq=floor(size(Var,1)/2);
117for 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
138end
139specmean=specmean/NbPos;
140
141%plot spectrum in log log
142hfig=findobj('Tag','fig_spectrum');
143if isempty(hfig)% create spectruim figure if it does not exist
144    hfig=figure;
145    set(hfig,'Tag','fig_spectrum');
146else
147    figure(hfig)
148end
149loglog(freq_max*(1:length(specmean))/length(specmean),specmean)
150set(gca,'YLim',[1.0000e-06*max(specmean) 1.1*max(specmean)])
151title (['power spectrum of ' Param.TransformInput.VariableName ])
152xlabel(['frequency (cycles per ' TimeUnit ')'])
153ylabel('spectral intensity')
154legend({'spectrum','cospectrum t t-1'})
155get(gca,'Unit')
156sum(specmean)
157if NbPos~=size(Var,2)
158    disp([ 'warning: ' num2str(size(Var,2)-NbPos) ' NaN sampled removed'])
159end
160if 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'])
163end
164grid on
165
166
Note: See TracBrowser for help on using the repository browser.