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

Last change on this file since 977 was 977, checked in by g7moreau, 4 years ago
  • Update Copyright 2008-2017 notice
File size: 8.3 KB
Line 
1% 'signal_spectrum': calculate and display power 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-2017, 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_spectrum(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 = {'nbre of points for the sliding window'};% 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        if DataOut.TransformInput.IndexRange(1,2)-DataOut.TransformInput.IndexRange(1,1)<DataOut.TransformInput.WindowLength
83            msgbox_uvmat('ERROR','WindowLength must be smaller than the total time index range')
84            return
85        end
86        huvmat=findobj(allchild(0),'Tag','uvmat');
87        UvData=get(huvmat,'UserData');
88        Data=UvData.PlotAxes;
89        YName=Data.ListVarName{1};
90        XName=Data.ListVarName{2};
91        yindex=DataOut.TransformInput.IndexRange(2,:);
92        y=Data.(YName)(yindex);
93        xindex=DataOut.TransformInput.IndexRange(3,:);
94        x=Data.(XName)(xindex);
95        haxes=findobj(huvmat,'Tag','PlotAxes');
96        axes(haxes);
97        hbounds=findobj(haxes,'Tag','Bounds');
98        if isempty(hbounds)
99        hbounds=rectangle('Position',[x(1) y(1) x(2)-x(1) y(2)-y(1)],'Tag','Bounds');
100        else
101            set(hbounds,'Position',[x(1) y(1) x(2)-x(1) y(2)-y(1)])
102        end
103    end
104    return
105end
106
107%% retrieve parameters
108DataOut=DataIn;
109WindowLength=Param.TransformInput.WindowLength;
110Shift=round(WindowLength/2);% shift between two windowsof analysis (half window length by default)
111
112%% get the variable to process
113if  ~isfield(DataIn,Param.TransformInput.VariableName)
114    return
115end
116Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse
117if isfield(Param.TransformInput,'IndexRange')
118    IndexRange=Param.TransformInput.IndexRange;
119    switch size(IndexRange,1)
120        case 3
121            Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2),IndexRange(3,1):IndexRange(3,2));
122        case 2
123            Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2));
124        case 1
125            Var=Var(IndexRange(1,1):IndexRange(1,2));
126    end
127end
128np=size(Var);%dimensions of Var
129if ~isvector(Var)
130    Var=reshape(Var,np(1),prod(np(2:end)));% reshape in a 2D matrix with time as first index
131end
132Var=Var-ones(np(1),1)*nanmean(Var,1); %substract mean value (excluding NaN)
133
134%% look for 'time' coordinate
135VarIndex=find(strcmp(Param.TransformInput.VariableName,DataIn.ListVarName));
136TimeDimName=DataIn.VarDimName{VarIndex}{1};
137TimeVarNameIndex=find(strcmp(TimeDimName,DataIn.ListVarName));
138if isempty(TimeVarNameIndex)
139    Time=1:np(1);
140    TimeUnit='vector index';
141else
142    Time=DataIn.(DataIn.ListVarName{TimeVarNameIndex});
143    TimeUnit=['Unit of ' TimeDimName];
144end
145% check time intervals
146diff_x=diff(Time);
147dx=min(diff_x);
148freq_max=1/(2*dx);
149check_interp=0;
150if diff_x>1.001*dx % non constant time interval
151    check_interp=1;
152end
153
154%% calculate the spectrum
155specmean=0;% mean spectrum initialisation
156cospecmean=0;
157NbNan=0;
158NbPos=0;
159for pos=1:size(Var,2)
160    sample=Var(:,pos);%extract sample to analyse
161    ind_bad=find(isnan(sample));
162    ind_good=find(~isnan(sample));
163    if numel(ind_good)>WindowLength
164        NbPos=NbPos+1;
165        if ~isempty(ind_bad)
166            sample=sample(ind_good); % keep only  non NaN data
167            NbNan=NbNan+numel(ind_bad);
168        end
169        %interpolate if needed
170        if ~isempty(ind_bad)||check_interp
171            sample=interp1(Time(ind_good),sample,(Time(1):dx:Time(end))); %interpolated func
172            sample(isnan(sample))=[];
173        end
174        spec=pwelch(sample,WindowLength);% calculate spectrum with Welch method
175        cospec=cpsd(sample,circshift(sample,[1 0]),WindowLength);% calculate the cospectrum with the sample shifted by 1 time unit
176        specmean=spec+specmean;
177        cospecmean=cospec+cospecmean;
178    end
179end
180specmean=specmean/NbPos;
181cospecmean=cospecmean/NbPos;
182
183%plot spectrum in log log
184hfig=findobj('Tag','fig_spectrum');
185if isempty(hfig)% create spectruim figure if it does not exist
186    hfig=figure;
187    set(hfig,'Tag','fig_spectrum');
188else
189    figure(hfig)
190end
191loglog(freq_max*(1:length(specmean))/length(specmean),specmean)
192hold on
193loglog(freq_max*(1:length(cospecmean))/length(cospecmean),cospecmean,'r')
194hold off
195title (['power spectrum of ' Param.TransformInput.VariableName ])
196xlabel(['frequency (cycles per ' TimeUnit ')'])
197ylabel('spectral intensity')
198legend({'spectrum','cospectrum t t-1'})
199get(gca,'Unit')
200sum(specmean)
201sum(cospecmean)
202if NbPos~=size(Var,2)
203    disp([ 'warning: ' num2str(size(Var,2)-NbPos) ' NaN sampled removed'])
204end
205if NbNan~=0
206    disp([ 'warning: ' num2str(NbNan) ' NaN values replaced by linear interpolation'])
207%text(0.9, 0.5,[ 'warning: ' num2str(NbNan) ' NaN values removed'])
208end
209grid on
210
211
Note: See TracBrowser for help on using the repository browser.