source: trunk/src/transform_field/signal_spectrum.m @ 1061

Last change on this file since 1061 was 1061, checked in by g7moreau, 5 years ago
  • Update copyright to year 2019
File size: 8.6 KB
RevLine 
[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)
[1010]3%  this function aplies the Welch method and calls the function  of the matlab signal processing toolbox
4%  it calculates also the cospectrum with the signal shifted by one
5%  index, to remove noise effects assumed uncorrelated from one iteration to
6%  the next.
[753]7%
8% OUTPUT:
9% DataOut: if DataIn.Action.RUN=0 (introducing parameters): Matlab structure containing the parameters
10%          else transformed field, here not modified (the function just produces a plot on an independent fig)
11%
12% INPUT:
13% DataIn: Matlab structure containing the input field from the GUI uvmat, DataIn.Action.RUN=0 to set input parameters.
14% Param: structure containing processing parameters, created when DataIn.Action.RUN=0 at the first use of the transform fct
15
[810]16%=======================================================================
[1061]17% Copyright 2008-2019, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
[810]18%   http://www.legi.grenoble-inp.fr
19%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
20%
21%     This file is part of the toolbox UVMAT.
22%
23%     UVMAT is free software; you can redistribute it and/or modify
24%     it under the terms of the GNU General Public License as published
25%     by the Free Software Foundation; either version 2 of the license,
26%     or (at your option) any later version.
27%
28%     UVMAT is distributed in the hope that it will be useful,
29%     but WITHOUT ANY WARRANTY; without even the implied warranty of
30%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
31%     GNU General Public License (see LICENSE.txt) for more details.
32%=======================================================================
33
[753]34function DataOut=signal_spectrum(DataIn,Param)
35
36%% request input parameters
37if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0)
38    VarNbDim=cellfun('length',DataIn.VarDimName);
39    [tild,rank]=sort(VarNbDim,2,'descend');% sort the list of input variables, putting the ones with higher dimensionality first
40    ListVarName=DataIn.ListVarName(rank);
41    VarDimName=DataIn.VarDimName(rank);
42    InitialValue=1;%default choice
43    if isfield(Param,'TransformInput') && isfield(Param.TransformInput,'VariableName')
44        val=find(strcmp(Param.TransformInput.VariableName,ListVarName));
45        if ~isempty(val);
46            InitialValue=val;
47        end
48    end
49    [s,OK] = listdlg('PromptString','Select the variable to process:',...
50        'SelectionMode','single','InitialValue',InitialValue,...
51        'ListString',ListVarName);
52    if OK==1
53        VarName=ListVarName{s};
54        DataOut.TransformInput.VariableName=VarName;
55        dlg_title = [mfilename ' calulates spectra along first dim ' VarDimName{s}{1}];% title of the input dialog fig
56        prompt = {'nbre of points for the sliding window'};% titles of the edit boxes
57        %default input:
58        def={'512'};% window length
59        np=size(DataIn.(VarName));
60        for idim=1:numel(np) % size restriction
61            if idim==1
62                prompt=[prompt;{['index range for spectral dim ' VarDimName{s}{idim}]}];% titles of the edit boxes
63            else
64            prompt=[prompt;{['index range for ' VarDimName{s}{idim}]}];% titles of the edit boxes
65            end
66            def=[def;{num2str([1 np(idim)])}];
67        end
68        if isfield(Param,'TransformInput')
69            if isfield(Param.TransformInput,'WindowLength')
70                def{1}=num2str(Param.TransformInput.WindowLength);
71            end
72            if isfield(Param.TransformInput,'IndexRange')
73                for ilist=1:min(numel(np),size(Param.TransformInput.IndexRange,1))
74                    def{ilist+1}=num2str(Param.TransformInput.IndexRange(ilist,:));
75                end
76            end
77        end
78        num_lines= 1;%numel(prompt);
79        % open the dialog fig
80        answer = inputdlg(prompt,dlg_title,num_lines,def);
81        DataOut.TransformInput.WindowLength=str2num(answer{1});
82        for ilist=1:numel(answer)-1
83            DataOut.TransformInput.IndexRange(ilist,1:2)=str2num(answer{ilist+1});
84        end
[944]85        if DataOut.TransformInput.IndexRange(1,2)-DataOut.TransformInput.IndexRange(1,1)<DataOut.TransformInput.WindowLength
86            msgbox_uvmat('ERROR','WindowLength must be smaller than the total time index range')
87            return
88        end
89        huvmat=findobj(allchild(0),'Tag','uvmat');
90        UvData=get(huvmat,'UserData');
91        Data=UvData.PlotAxes;
92        YName=Data.ListVarName{1};
93        XName=Data.ListVarName{2};
94        yindex=DataOut.TransformInput.IndexRange(2,:);
95        y=Data.(YName)(yindex);
[1010]96        if size(DataOut.TransformInput.IndexRange,1)>=3
[944]97        xindex=DataOut.TransformInput.IndexRange(3,:);
[1010]98        else
99            xindex=DataOut.TransformInput.IndexRange(2,:);
100        end
[944]101        x=Data.(XName)(xindex);
102        haxes=findobj(huvmat,'Tag','PlotAxes');
103        axes(haxes);
104        hbounds=findobj(haxes,'Tag','Bounds');
105        if isempty(hbounds)
106        hbounds=rectangle('Position',[x(1) y(1) x(2)-x(1) y(2)-y(1)],'Tag','Bounds');
107        else
108            set(hbounds,'Position',[x(1) y(1) x(2)-x(1) y(2)-y(1)])
109        end
[753]110    end
111    return
112end
113
114%% retrieve parameters
115DataOut=DataIn;
116WindowLength=Param.TransformInput.WindowLength;
117Shift=round(WindowLength/2);% shift between two windowsof analysis (half window length by default)
118
119%% get the variable to process
[944]120if  ~isfield(DataIn,Param.TransformInput.VariableName)
121    return
122end
[753]123Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse
[759]124if isfield(Param.TransformInput,'IndexRange')
125    IndexRange=Param.TransformInput.IndexRange;
126    switch size(IndexRange,1)
127        case 3
128            Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2),IndexRange(3,1):IndexRange(3,2));
129        case 2
130            Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2));
131        case 1
132            Var=Var(IndexRange(1,1):IndexRange(1,2));
133    end
134end
[753]135np=size(Var);%dimensions of Var
136if ~isvector(Var)
137    Var=reshape(Var,np(1),prod(np(2:end)));% reshape in a 2D matrix with time as first index
138end
[759]139Var=Var-ones(np(1),1)*nanmean(Var,1); %substract mean value (excluding NaN)
[753]140
141%% look for 'time' coordinate
142VarIndex=find(strcmp(Param.TransformInput.VariableName,DataIn.ListVarName));
143TimeDimName=DataIn.VarDimName{VarIndex}{1};
144TimeVarNameIndex=find(strcmp(TimeDimName,DataIn.ListVarName));
145if isempty(TimeVarNameIndex)
146    Time=1:np(1);
147    TimeUnit='vector index';
148else
149    Time=DataIn.(DataIn.ListVarName{TimeVarNameIndex});
150    TimeUnit=['Unit of ' TimeDimName];
151end
152% check time intervals
153diff_x=diff(Time);
154dx=min(diff_x);
155freq_max=1/(2*dx);
156check_interp=0;
157if diff_x>1.001*dx % non constant time interval
158    check_interp=1;
159end
160
[944]161%% calculate the spectrum
[753]162specmean=0;% mean spectrum initialisation
163cospecmean=0;
164NbNan=0;
165NbPos=0;
166for pos=1:size(Var,2)
167    sample=Var(:,pos);%extract sample to analyse
168    ind_bad=find(isnan(sample));
169    ind_good=find(~isnan(sample));
170    if numel(ind_good)>WindowLength
171        NbPos=NbPos+1;
172        if ~isempty(ind_bad)
173            sample=sample(ind_good); % keep only  non NaN data
174            NbNan=NbNan+numel(ind_bad);
175        end
176        %interpolate if needed
177        if ~isempty(ind_bad)||check_interp
178            sample=interp1(Time(ind_good),sample,(Time(1):dx:Time(end))); %interpolated func
[759]179            sample(isnan(sample))=[];
[753]180        end
181        spec=pwelch(sample,WindowLength);% calculate spectrum with Welch method
182        cospec=cpsd(sample,circshift(sample,[1 0]),WindowLength);% calculate the cospectrum with the sample shifted by 1 time unit
183        specmean=spec+specmean;
184        cospecmean=cospec+cospecmean;
185    end
186end
187specmean=specmean/NbPos;
188cospecmean=cospecmean/NbPos;
189
190%plot spectrum in log log
191hfig=findobj('Tag','fig_spectrum');
192if isempty(hfig)% create spectruim figure if it does not exist
193    hfig=figure;
194    set(hfig,'Tag','fig_spectrum');
195else
196    figure(hfig)
197end
198loglog(freq_max*(1:length(specmean))/length(specmean),specmean)
199hold on
200loglog(freq_max*(1:length(cospecmean))/length(cospecmean),cospecmean,'r')
201hold off
202title (['power spectrum of ' Param.TransformInput.VariableName ])
203xlabel(['frequency (cycles per ' TimeUnit ')'])
204ylabel('spectral intensity')
205legend({'spectrum','cospectrum t t-1'})
206get(gca,'Unit')
[759]207sum(specmean)
208sum(cospecmean)
[753]209if NbPos~=size(Var,2)
210    disp([ 'warning: ' num2str(size(Var,2)-NbPos) ' NaN sampled removed'])
211end
212if NbNan~=0
213    disp([ 'warning: ' num2str(NbNan) ' NaN values replaced by linear interpolation'])
214%text(0.9, 0.5,[ 'warning: ' num2str(NbNan) ' NaN values removed'])
215end
216grid on
217
218
Note: See TracBrowser for help on using the repository browser.