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

Last change on this file since 933 was 924, checked in by g7moreau, 9 years ago
  • Update Copyright Copyright 2008-2016, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
File size: 7.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-2016, 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    end
83    return
84end
85
86%% retrieve parameters
87DataOut=DataIn;
88WindowLength=Param.TransformInput.WindowLength;
89Shift=round(WindowLength/2);% shift between two windowsof analysis (half window length by default)
90
91%% get the variable to process
92Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse
93if isfield(Param.TransformInput,'IndexRange')
94    IndexRange=Param.TransformInput.IndexRange;
95    switch size(IndexRange,1)
96        case 3
97            Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2),IndexRange(3,1):IndexRange(3,2));
98        case 2
99            Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2));
100        case 1
101            Var=Var(IndexRange(1,1):IndexRange(1,2));
102    end
103end
104np=size(Var);%dimensions of Var
105if ~isvector(Var)
106    Var=reshape(Var,np(1),prod(np(2:end)));% reshape in a 2D matrix with time as first index
107end
108Var=Var-ones(np(1),1)*nanmean(Var,1); %substract mean value (excluding NaN)
109
110%% look for 'time' coordinate
111VarIndex=find(strcmp(Param.TransformInput.VariableName,DataIn.ListVarName));
112TimeDimName=DataIn.VarDimName{VarIndex}{1};
113TimeVarNameIndex=find(strcmp(TimeDimName,DataIn.ListVarName));
114if isempty(TimeVarNameIndex)
115    Time=1:np(1);
116    TimeUnit='vector index';
117else
118    Time=DataIn.(DataIn.ListVarName{TimeVarNameIndex});
119    TimeUnit=['Unit of ' TimeDimName];
120end
121% check time intervals
122diff_x=diff(Time);
123dx=min(diff_x);
124freq_max=1/(2*dx);
125check_interp=0;
126if diff_x>1.001*dx % non constant time interval
127    check_interp=1;
128end
129
130%% claculate the spectrum
131specmean=0;% mean spectrum initialisation
132cospecmean=0;
133NbNan=0;
134NbPos=0;
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        spec=pwelch(sample,WindowLength);% calculate spectrum with Welch method
151        cospec=cpsd(sample,circshift(sample,[1 0]),WindowLength);% calculate the cospectrum with the sample shifted by 1 time unit
152        specmean=spec+specmean;
153        cospecmean=cospec+cospecmean;
154    end
155end
156specmean=specmean/NbPos;
157cospecmean=cospecmean/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)
168hold on
169loglog(freq_max*(1:length(cospecmean))/length(cospecmean),cospecmean,'r')
170hold off
171title (['power spectrum of ' Param.TransformInput.VariableName ])
172xlabel(['frequency (cycles per ' TimeUnit ')'])
173ylabel('spectral intensity')
174legend({'spectrum','cospectrum t t-1'})
175get(gca,'Unit')
176sum(specmean)
177sum(cospecmean)
178if NbPos~=size(Var,2)
179    disp([ 'warning: ' num2str(size(Var,2)-NbPos) ' NaN sampled removed'])
180end
181if NbNan~=0
182    disp([ 'warning: ' num2str(NbNan) ' NaN values replaced by linear interpolation'])
183%text(0.9, 0.5,[ 'warning: ' num2str(NbNan) ' NaN values removed'])
184end
185grid on
186
187
Note: See TracBrowser for help on using the repository browser.