Index: /trunk/src/transform_field/ima_filter.m =================================================================== --- /trunk/src/transform_field/ima_filter.m (revision 752) +++ /trunk/src/transform_field/ima_filter.m (revision 753) @@ -12,6 +12,8 @@ %------------------------------------- -function DataOut=ima_test(DataIn,Param) -if isequal(DataIn,'*') +function DataOut=ima_filter(DataIn,Param) + +%% request input parameters +if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0) prompt = {'npx';'npy'}; dlg_title = 'get the filter size in x and y'; Index: /trunk/src/transform_field/phys.m =================================================================== --- /trunk/src/transform_field/phys.m (revision 752) +++ /trunk/src/transform_field/phys.m (revision 753) @@ -28,5 +28,5 @@ DataOut=[]; DataOut_1=[]; %default second output field -if strcmp(DataIn,'*') +if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0) if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit') DataOut.CoordUnit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units Index: /trunk/src/transform_field/signal_bandpass_filter.m =================================================================== --- /trunk/src/transform_field/signal_bandpass_filter.m (revision 753) +++ /trunk/src/transform_field/signal_bandpass_filter.m (revision 753) @@ -0,0 +1,30 @@ +% 'signal_band_filter': example of pass-band filter (using the signal toolbox fct fdesign) +% frequency and bandwidth need to be modified in the function +% +% OUTPUT: +% DataOut: Matlab structure representing the output (filtered) field +% +%INPUT: +% DataIn: Matlab structure representing the output field +% + +function DataOut=signal_bandpass_filter(DataIn) +%% set GUI config +DataOut=[]; +if strcmp(DataIn,'*') + DataOut.InputFieldType='1D'; + return +end +%%%%%%%%%%%%%%%%%%%%%%%%% +frequency=2*5/54.5;% frequency at which the signal needs to be filetered + %d=fdesign.bandpass(0.8*frequency, 0.95*frequency, 1.05*frequency, 1.2*frequency, 60, 1, 60); + d=fdesign.bandpass(0.7*frequency, 0.9*frequency, 1.1*frequency, 1.3*frequency, 60, 1, 60); + Hd=design(d);% command fvtool(Hd) to visualize the filter frequency response + DataOut=DataIn; + if isfield(DataIn,'U') + DataOut.U = filter(Hd,DataIn.U); + end + if isfield(DataIn,'V') + DataOut.V = filter(Hd,DataIn.V); + end + Index: unk/src/transform_field/signal_filter.m =================================================================== --- /trunk/src/transform_field/signal_filter.m (revision 752) +++ (revision ) @@ -1,31 +1,0 @@ -% 'FFT': calculate and display 2D spectrum of the input scalar -% GUI_input=FFT(hget_field) -% -% OUTPUT: -% GUI_input: option for display in the GUI get_field -% -%INPUT: -% hget_field: handles of the GUI get_field -% - -function DataOut=signal_filter(DataIn) -%% set GUI config -DataOut=[]; -if strcmp(DataIn,'*') - DataOut.InputFieldType='1D'; - return -end -%%%%%%%%%%%%%%%%%%%%%%%%% -frequency=2*5/54.5; -frequency=frequency; - %d=fdesign.bandpass(0.8*frequency, 0.95*frequency, 1.05*frequency, 1.2*frequency, 60, 1, 60); - d=fdesign.bandpass(0.7*frequency, 0.9*frequency, 1.1*frequency, 1.3*frequency, 60, 1, 60); - Hd=design(d);% command fvtool(Hd) to visualize the filter frequency response - DataOut=DataIn; - if isfield(DataIn,'U') - DataOut.U = filter(Hd,DataIn.U); - end - if isfield(DataIn,'V') - DataOut.V = filter(Hd,DataIn.V); - end - Index: /trunk/src/transform_field/signal_spectrum.m =================================================================== --- /trunk/src/transform_field/signal_spectrum.m (revision 753) +++ /trunk/src/transform_field/signal_spectrum.m (revision 753) @@ -0,0 +1,155 @@ +% '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_spectrum(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 = {'nbre of points for the sliding window'};% 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; +Shift=round(WindowLength/2);% shift between two windowsof analysis (half window length by default) + +%% get the variable to process +Var= DataIn.(Param.TransformInput.VariableName);%variable to analyse +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; +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 + end + spec=pwelch(sample,WindowLength);% calculate spectrum with Welch method + cospec=cpsd(sample,circshift(sample,[1 0]),WindowLength);% calculate the cospectrum with the sample shifted by 1 time unit + specmean=spec+specmean; + cospecmean=cospec+cospecmean; + end +end +specmean=specmean/NbPos; +cospecmean=cospecmean/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) +hold on +loglog(freq_max*(1:length(cospecmean))/length(cospecmean),cospecmean,'r') +hold off +title (['power spectrum of ' Param.TransformInput.VariableName ]) +xlabel(['frequency (cycles per ' TimeUnit ')']) +ylabel('spectral intensity') +legend({'spectrum','cospectrum t t-1'}) +get(gca,'Unit') +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 + +