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

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

signal spectrum analysis and various improvements in transform

File size: 6.0 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_spectrum(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 = {'nbre of points for the sliding window'};% 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;
71Shift=round(WindowLength/2);% shift between two windowsof analysis (half window length by default)
72
73%% get the variable to process
74Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse
75np=size(Var);%dimensions of Var
76if ~isvector(Var)
77    Var=reshape(Var,np(1),prod(np(2:end)));% reshape in a 2D matrix with time as first index
78end
79Var=Var-ones(np(1),1)*nanmean(Var,1); %substract mean value (excluding NaN)
80
81%% look for 'time' coordinate
82VarIndex=find(strcmp(Param.TransformInput.VariableName,DataIn.ListVarName));
83TimeDimName=DataIn.VarDimName{VarIndex}{1};
84TimeVarNameIndex=find(strcmp(TimeDimName,DataIn.ListVarName));
85if isempty(TimeVarNameIndex)
86    Time=1:np(1);
87    TimeUnit='vector index';
88else
89    Time=DataIn.(DataIn.ListVarName{TimeVarNameIndex});
90    TimeUnit=['Unit of ' TimeDimName];
91end
92% check time intervals
93diff_x=diff(Time);
94dx=min(diff_x);
95freq_max=1/(2*dx);
96check_interp=0;
97if diff_x>1.001*dx % non constant time interval
98    check_interp=1;
99end
100
101%% claculate the spectrum
102specmean=0;% mean spectrum initialisation
103cospecmean=0;
104NbNan=0;
105NbPos=0;
106for pos=1:size(Var,2)
107    sample=Var(:,pos);%extract sample to analyse
108    ind_bad=find(isnan(sample));
109    ind_good=find(~isnan(sample));
110    if numel(ind_good)>WindowLength
111        NbPos=NbPos+1;
112        if ~isempty(ind_bad)
113            sample=sample(ind_good); % keep only  non NaN data
114            NbNan=NbNan+numel(ind_bad);
115        end
116        %interpolate if needed
117        if ~isempty(ind_bad)||check_interp
118            sample=interp1(Time(ind_good),sample,(Time(1):dx:Time(end))); %interpolated func
119        end
120        spec=pwelch(sample,WindowLength);% calculate spectrum with Welch method
121        cospec=cpsd(sample,circshift(sample,[1 0]),WindowLength);% calculate the cospectrum with the sample shifted by 1 time unit
122        specmean=spec+specmean;
123        cospecmean=cospec+cospecmean;
124    end
125end
126specmean=specmean/NbPos;
127cospecmean=cospecmean/NbPos;
128
129%plot spectrum in log log
130hfig=findobj('Tag','fig_spectrum');
131if isempty(hfig)% create spectruim figure if it does not exist
132    hfig=figure;
133    set(hfig,'Tag','fig_spectrum');
134else
135    figure(hfig)
136end
137loglog(freq_max*(1:length(specmean))/length(specmean),specmean)
138hold on
139loglog(freq_max*(1:length(cospecmean))/length(cospecmean),cospecmean,'r')
140hold off
141title (['power spectrum of ' Param.TransformInput.VariableName ])
142xlabel(['frequency (cycles per ' TimeUnit ')'])
143ylabel('spectral intensity')
144legend({'spectrum','cospectrum t t-1'})
145get(gca,'Unit')
146if NbPos~=size(Var,2)
147    disp([ 'warning: ' num2str(size(Var,2)-NbPos) ' NaN sampled removed'])
148end
149if NbNan~=0
150    disp([ 'warning: ' num2str(NbNan) ' NaN values replaced by linear interpolation'])
151%text(0.9, 0.5,[ 'warning: ' num2str(NbNan) ' NaN values removed'])
152end
153grid on
154
155
Note: See TracBrowser for help on using the repository browser.