source: trunk/src/get_field/FFT.m @ 146

Last change on this file since 146 was 93, checked in by sommeria, 14 years ago

FFT: improved to deal with NaN data
merge_proj: corrected to reproduce dt if unique
uvmat.fig: tooltip corrected
plot_field: bug for isocontour corrected,
im_filter: cleaning
phys_polar: spatial derivative included (still to check)
set_obeject.fig: minor correction
struct2nc: comments improved
uvmat: button NB implemented
read_civxdata: error message improved

File size: 2.9 KB
Line 
1function GUI_input=FFT(hget_field)
2global spec x_vec
3%requests for the visibility of input windows in the GUI series  (activated directly by the selection in the menu ACTION)
4if ~exist('hget_field','var')
5    GUI_input={'check_1Dplot'};
6    return %exit the function
7end
8GUI_input=[];
9%initiation
10hhget_field=guidata(hget_field);
11% testinterp=0;
12abscissa_list=get(hhget_field.abscissa,'String');
13val=get(hhget_field.abscissa,'Value');
14val=val(1);
15abscissa_name=abscissa_list{val};
16ordinate_list=get(hhget_field.ordinate,'String');
17val=get(hhget_field.ordinate,'Value');
18val=val(1); %take only the first variable in the list
19
20%ordinate_name=Field.ListVarName{val};
21ordinate_name=ordinate_list{val};
22
23[Field,errormsg]=read_get_field(hget_field);
24if ~isempty(errormsg)
25    msgbox_uvmat('ERROR',['error in get_field/FFT input:' errormsg])
26    return
27end
28
29% get variable
30eval(['Var= Field.' ordinate_name ';']);
31np=size(Var)
32np_freq=floor(np(1)/2);
33dx=1;%default
34dfreq=1/np(1);%default frequency interval (abscissa= array index)
35sum_data=sum(Var,2);
36if ~isequal(abscissa_name,'')
37    eval(['Coord_x= Field.' abscissa_name ';']);
38    ind_select=find(~isinf(Coord_x)&~isnan(sum_data));%detect infinite values
39    Coord_x=Coord_x(ind_select);
40    Var=Var(ind_select,:);
41    diff_x=diff(Coord_x);
42    dx=min(diff_x);
43    %interpolate on a regular abscissa interval if needed
44    if (max(diff_x)-dx)> 0.001*dx || numel(ind_select)<np(1)
45        xequ=Coord_x(1):dx:Coord_x(end);%equal time spacingdx=
46        Var=interp1(Coord_x,Var,xequ); %interpolated func
47        np=size(Var);
48    end
49 %   funcinterp=interp1(time,func,timeq); %interpolated func
50    dfreq=1/(Coord_x(end)-Coord_x(1));%frequency interval
51end
52freq_max=1/(2*dx);
53Var=Var-ones(np(1),1)*mean(Var,1); %substract mean value
54fourier=fft(Var);%take fft (complex)
55spec=abs(fourier).*abs(fourier);% take square of the modulus
56spec=spec(1:np_freq,:);%keep only the first half (the other is symmetric)
57
58%plot
59list_fig=get(hhget_field.list_fig,'String');
60val=get(hhget_field.list_fig,'Value');
61hfig=str2num(list_fig{val})% chosen figure number from tyhe GUI
62if isempty(hfig)
63    hfig=figure;
64else
65    figure(hfig);
66end
67haxes=findobj(hfig,'Type','axes');
68if ~isempty(haxes)
69    axes(haxes)
70end
71x_vec=linspace(dfreq,freq_max,np_freq);
72plot(x_vec',spec)
73xlabel('frequency (Hz)')
74ylabel('spectral intensity')
75grid on
76
77%
78%
79% np=length(funcinterp);
80% funcinterp=funcinterp-sum(funcinterp)/np; %substract mean
81% fourier=fft(funcinterp);%take fft (complex)
82% spec=abs(fourier).*abs(fourier);% take sqare of the modulus
83% spec=spec([1:floor(np/2)]);%keep only the first half (the other is symmetric)
84% eval(['Field.' varname '=spec;'])
85% Field
86% % dfreq=1/(time(end)-time(1));%frequency interval
87% % freq=[0:dfreq:(floor(np/2)-1)*dfreq];
88% % figure(1)
89% % hold on
90% % plot(freq,spec)
91% % xlabel('frequency (Hz)')
92% % ylabel('spectral intensity')
93% % title(['spectrum of' fields]);
94% % grid on
95
Note: See TracBrowser for help on using the repository browser.