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

Last change on this file since 1040 was 1027, checked in by g7moreau, 6 years ago
  • Update Copyright 2017 -> 2018
File size: 7.1 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
13%=======================================================================
14% Copyright 2008-2018, LEGI UMR 5519 / CNRS UGA 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
31function DataOut=signal_FFTMean(DataIn,Param)
32
33%% request input parameters
34if 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
84end
85
86%% retrieve parameters
87DataOut=DataIn;
88WindowLength=Param.TransformInput.WindowLength;
89
90%% get the variable to process
91Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse
92if 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
102end
103np=size(Var);%dimensions of Var
104if ~isvector(Var)
105    Var=reshape(Var,np(1),prod(np(2:end)));% reshape in a 2D matrix with time as first index
106end
107Var=Var-ones(np(1),1)*nanmean(Var,1); %substract mean value (excluding NaN)
108
109%% look for 'time' coordinate
110VarIndex=find(strcmp(Param.TransformInput.VariableName,DataIn.ListVarName));
111TimeDimName=DataIn.VarDimName{VarIndex}{1};
112TimeVarNameIndex=find(strcmp(TimeDimName,DataIn.ListVarName));
113if isempty(TimeVarNameIndex)
114    Time=1:np(1);
115    TimeUnit='vector index';
116else
117    Time=DataIn.(DataIn.ListVarName{TimeVarNameIndex});
118    TimeUnit=['Unit of ' TimeDimName];
119end
120% check time intervals
121diff_x=diff(Time);
122dx=min(diff_x);
123freq_max=1/(2*dx);
124check_interp=0;
125if diff_x>1.001*dx % non constant time interval
126    check_interp=1;
127end
128
129%% claculate the spectrum
130specmean=0;% mean spectrum initialisation
131cospecmean=0;
132NbNan=0;
133NbPos=0;
134np_freq=floor(size(Var,1)/2);
135for 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
156end
157specmean=specmean/NbPos;
158
159%plot spectrum in log log
160hfig=findobj('Tag','fig_spectrum');
161if isempty(hfig)% create spectruim figure if it does not exist
162    hfig=figure;
163    set(hfig,'Tag','fig_spectrum');
164else
165    figure(hfig)
166end
167loglog(freq_max*(1:length(specmean))/length(specmean),specmean)
168set(gca,'YLim',[1.0000e-06*max(specmean) 1.1*max(specmean)])
169title (['power spectrum of ' Param.TransformInput.VariableName ])
170xlabel(['frequency (cycles per ' TimeUnit ')'])
171ylabel('spectral intensity')
172legend({'spectrum','cospectrum t t-1'})
173get(gca,'Unit')
174sum(specmean)
175if NbPos~=size(Var,2)
176    disp([ 'warning: ' num2str(size(Var,2)-NbPos) ' NaN sampled removed'])
177end
178if 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'])
181end
182grid on
183
184
Note: See TracBrowser for help on using the repository browser.