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