Index: trunk/src/transform_field/signal_FFTMean.m =================================================================== --- trunk/src/transform_field/signal_FFTMean.m (revision 766) +++ trunk/src/transform_field/signal_FFTMean.m (revision 766) @@ -0,0 +1,166 @@ +% 'signal_spectrum': calculate and display spectrum of the current field +% operate on a 1D signal or the first dimension of a higher dimensional matrix (then average over other dimensions) +% this function aplies the Welch method and call the function of the matlab signal processing toolbox +% +% OUTPUT: +% DataOut: if DataIn.Action.RUN=0 (introducing parameters): Matlab structure containing the parameters +% else transformed field, here not modified (the function just produces a plot on an independent fig) +% +% INPUT: +% DataIn: Matlab structure containing the input field from the GUI uvmat, DataIn.Action.RUN=0 to set input parameters. +% Param: structure containing processing parameters, created when DataIn.Action.RUN=0 at the first use of the transform fct + +function DataOut=signal_FFTMean(DataIn,Param) + +%% request input parameters +if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0) + VarNbDim=cellfun('length',DataIn.VarDimName); + [tild,rank]=sort(VarNbDim,2,'descend');% sort the list of input variables, putting the ones with higher dimensionality first + ListVarName=DataIn.ListVarName(rank); + VarDimName=DataIn.VarDimName(rank); + InitialValue=1;%default choice + if isfield(Param,'TransformInput') && isfield(Param.TransformInput,'VariableName') + val=find(strcmp(Param.TransformInput.VariableName,ListVarName)); + if ~isempty(val); + InitialValue=val; + end + end + [s,OK] = listdlg('PromptString','Select the variable to process:',... + 'SelectionMode','single','InitialValue',InitialValue,... + 'ListString',ListVarName); + if OK==1 + VarName=ListVarName{s}; + DataOut.TransformInput.VariableName=VarName; + dlg_title = [mfilename ' calulates spectra along first dim ' VarDimName{s}{1}];% title of the input dialog fig + prompt = {'not used'};% titles of the edit boxes + %default input: + def={'512'};% window length + np=size(DataIn.(VarName)); + for idim=1:numel(np) % size restriction + if idim==1 + prompt=[prompt;{['index range for spectral dim ' VarDimName{s}{idim}]}];% titles of the edit boxes + else + prompt=[prompt;{['index range for ' VarDimName{s}{idim}]}];% titles of the edit boxes + end + def=[def;{num2str([1 np(idim)])}]; + end + if isfield(Param,'TransformInput') + if isfield(Param.TransformInput,'WindowLength') + def{1}=num2str(Param.TransformInput.WindowLength); + end + if isfield(Param.TransformInput,'IndexRange') + for ilist=1:min(numel(np),size(Param.TransformInput.IndexRange,1)) + def{ilist+1}=num2str(Param.TransformInput.IndexRange(ilist,:)); + end + end + end + num_lines= 1;%numel(prompt); + % open the dialog fig + answer = inputdlg(prompt,dlg_title,num_lines,def); + DataOut.TransformInput.WindowLength=str2num(answer{1}); + for ilist=1:numel(answer)-1 + DataOut.TransformInput.IndexRange(ilist,1:2)=str2num(answer{ilist+1}); + end + end + return +end + +%% retrieve parameters +DataOut=DataIn; +WindowLength=Param.TransformInput.WindowLength; + +%% get the variable to process +Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse +if isfield(Param.TransformInput,'IndexRange') + IndexRange=Param.TransformInput.IndexRange; + switch size(IndexRange,1) + case 3 + Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2),IndexRange(3,1):IndexRange(3,2)); + case 2 + Var=Var(IndexRange(1,1):IndexRange(1,2),IndexRange(2,1):IndexRange(2,2)); + case 1 + Var=Var(IndexRange(1,1):IndexRange(1,2)); + end +end +np=size(Var);%dimensions of Var +if ~isvector(Var) + Var=reshape(Var,np(1),prod(np(2:end)));% reshape in a 2D matrix with time as first index +end +Var=Var-ones(np(1),1)*nanmean(Var,1); %substract mean value (excluding NaN) + +%% look for 'time' coordinate +VarIndex=find(strcmp(Param.TransformInput.VariableName,DataIn.ListVarName)); +TimeDimName=DataIn.VarDimName{VarIndex}{1}; +TimeVarNameIndex=find(strcmp(TimeDimName,DataIn.ListVarName)); +if isempty(TimeVarNameIndex) + Time=1:np(1); + TimeUnit='vector index'; +else + Time=DataIn.(DataIn.ListVarName{TimeVarNameIndex}); + TimeUnit=['Unit of ' TimeDimName]; +end +% check time intervals +diff_x=diff(Time); +dx=min(diff_x); +freq_max=1/(2*dx); +check_interp=0; +if diff_x>1.001*dx % non constant time interval + check_interp=1; +end + +%% claculate the spectrum +specmean=0;% mean spectrum initialisation +cospecmean=0; +NbNan=0; +NbPos=0; +np_freq=floor(size(Var,1)/2); +for pos=1:size(Var,2) + sample=Var(:,pos);%extract sample to analyse + ind_bad=find(isnan(sample)); + ind_good=find(~isnan(sample)); +% if numel(ind_good)>WindowLength + NbPos=NbPos+1; + if ~isempty(ind_bad) + sample=sample(ind_good); % keep only non NaN data + NbNan=NbNan+numel(ind_bad); + end + %interpolate if needed + if ~isempty(ind_bad)||check_interp + sample=interp1(Time(ind_good),sample,(Time(1):dx:Time(end))); %interpolated func + sample(isnan(sample))=[]; + end + + fourier=fft(sample);%take fft (complex) + spec=abs(fourier).*abs(fourier);% take square of the modulus + spec=spec(1:np_freq,:);%keep only the first half (the other is symmetric) + specmean=spec+specmean; +% end +end +specmean=specmean/NbPos; + +%plot spectrum in log log +hfig=findobj('Tag','fig_spectrum'); +if isempty(hfig)% create spectruim figure if it does not exist + hfig=figure; + set(hfig,'Tag','fig_spectrum'); +else + figure(hfig) +end +loglog(freq_max*(1:length(specmean))/length(specmean),specmean) +set(gca,'YLim',[1.0000e-06*max(specmean) 1.1*max(specmean)]) +title (['power spectrum of ' Param.TransformInput.VariableName ]) +xlabel(['frequency (cycles per ' TimeUnit ')']) +ylabel('spectral intensity') +legend({'spectrum','cospectrum t t-1'}) +get(gca,'Unit') +sum(specmean) +if NbPos~=size(Var,2) + disp([ 'warning: ' num2str(size(Var,2)-NbPos) ' NaN sampled removed']) +end +if NbNan~=0 + disp([ 'warning: ' num2str(NbNan) ' NaN values replaced by linear interpolation']) +%text(0.9, 0.5,[ 'warning: ' num2str(NbNan) ' NaN values removed']) +end +grid on + +