Two-dimensional spectral analysis Contents Make a synthesized 2-D signal Lets take the powerspectrum in x and look at it... Taking the 2-D Power Spectrum.... A more interesting process. Do the 2-D spectral estimate... Make a synthesized 2-D signal Here, let $ and $ be the space the data occupies, and $ be the signal. xx = 1:1000; yy = 1:2000; =meshgrid(xx,yy); kx0 = 0.01; ky0 = 0.004; w = cos(x*2*pi*kx0).*cos(y*2*pi*ky0); figure(1);clf imagesc(xx,yy,w); set(gca, 'ydir' , 'no' ); xlabel( 'x ' ); ylabel( 'y ' ); Lets take the powerspectrum in x and look at it... % first k_x = /\\deltax M = length(xx); N = length(yy); dx = 1; k_x = (0:M-1)/dx/M; dy = 1; k_y = (0:N-1)/dx/N; X = fft(w\')\'; Gx = (1/M/dx)*X(:,2:M/2).* conj(X(:,2:M/2)); k_x = k_x(2:M/2); The spectrum of any individual row is just the spectrum of w at that value of y: For this example, it is just a spike at $. figure(2);clf loglog(k_x,Gx(10,:)) xlabel( 'k_x ' ); ylabel( 'G_x(k_x: y=10 m)' ); The value of this peak varies with y as $. in = find(k_xkx0); in = in(1)-1; figure(3); clf; plot(y,Gx(:,in)); xlabel( 'y ' ); ylabel( 'G_x(k_x=k_{x0},y)' ); or visualized in 2-D.... figure(4); pcolor(k_x,yy,real(Gx));shading flat ; set(gca, 'xscale' , 'log' ); xlabel( 'k_x ' ); ylabel( 'y ' ); Taking the 2-D Power Spectrum.... XY = fft(fft(w\')\'); k_x = (0:M-1)/dx/M; k_y = (0:N-1)/dx/N; GXY = (1/M/dx)*(1/N/dy)*XY(2:N/2,2:M/2).* conj(XY(2:N/2,2:M/2)); k_y = k_y(2:N/2); k_x = k_x(2:M/2); We see the result is a single peak at $ figure(5);clf subplot(1,2,1); pcolor(k_x,k_y,real(GXY));shading flat ; set(gca, 'xscale' , 'log' , 'yscale' , 'log' ) xlabel( 'k_x ' ); ylabel( 'k_y ' ); subplot(1,2,2); surface(k_x,k_y,real(GXY));shading interp ; set(gca, 'xscale' , 'log' , 'yscale' , 'log' ) xlabel( 'k_x ' ); ylabel( 'k_y ' ); view(30,30); w0 = w; A more interesting process. Consider a 2-D white noise process, and then integrate separately in both x and y to get red processes in both dimensions.. Then ad the noise back on... wp = randn(N,M); w = cumsum(cumsum(wp)\')\'; w = w+wp*50+1000*w0; figure(10);clf subplot(2,1,1); imagesc(xx,yy,wp); set(gca, 'ydir' , 'nor' ); title( 'w''' ); ylabel( 'y ' ); subplot(2,1,2); imagesc(xx,yy,w); title( 'w' ); set(gca, 'ydir' , 'nor' ); xlabel( 'x ' ); Do the 2-D spectral estimate... XY = fft(fft(w\')\'); k_x = (0:M-1)/dx/M; k_y = (0:N-1)/dx/N; GXY = (1/M/dx)*(1/N/dy)*XY(2:N/2,2:M/2).* conj(XY(2:N/2,2:M/2)); k_y = k_y(2:N/2); k_x = k_x(2:M/2); Note how the spectrum has large amlitude at low wavenumbers, and then hits a noise floor at high wavenumbers. figure(15);clf subplot(1,2,1); pcolor(k_x,k_y,log10(real(GXY)));shading flat ; set(gca, 'xscale' , 'log' , 'yscale' , 'log' ) xlabel( 'k_x ' ); ylabel( 'k_y ' ); subplot(1,2,2); surface(k_x,k_y,log10(real(GXY)));shading interp ; set(gca, 'xscale' , 'log' , 'yscale' , 'log' ) xlabel( 'k_x ' ); ylabel( 'k_y ' ); view(40,20); Published with MATLAB® 7.7 原文请访问： http://valdez.seos.uvic.ca/~jklymak/Phy580/html/Fft2d.html
Problems w/ power spectra? Contents Getting started Aliasing Finite record length Spectral Leakage: Next Getting started load SynthSig ; Aliasing Here we need to generate another peak. om_peak= 4.4*2*pi; ph = rand(1,1)*2*pi; sig.x_peak2 = 2*10.*cos(om_peak*sig.t+ph)/10; =rawPowerSpec(sig.x_peak2+sig.x_peak+sig.x_w,dt); % Now, suppose that we can only sample at 0.25 Hz (i.e. 1/4 of the time) xx = sig.x_peak2+sig.x_peak+sig.x_w; x=xx(1:4:end); t = sig.t(1:4:end); =rawPowerSpec(x,median(diff(t))); jmkfigure(100,2,0.6);clf subplot(2,1,1); plot(x(1:100)); hold on ; plot(1:4:100,x(1:4:100), 'b.' , 'markersi' ,10); xlabel( 't ' ); ylabel( 'x ' ); subplot(2,1,2); loglog(f,p,fn,pn) axis tight ; yl = get(gca, 'ylim' ); xlabel( 'f ' ); ylabel( '\\phi_x ' ); axis tight print -dpng -r250 doc/Aliasing here we see the result of aliasing , where the peak has wrapped around to a harmonic in the resolved spectrum. To stop this we should smooth the new signal first: =butter(2,0.3); x = filtfilt(b,a,xx); subplot(2,1,1); plot(1:4:100,x(1:4:100), 'r.' , 'markersize' ,10) plot(x(1:100), 'g' ) xlabel( 't ' ); ylabel( 'x ' ); =rawPowerSpec(x,median(diff(sig.t))); x = x(1:4:end); t = sig.t(1:4:end); =rawPowerSpec(x,median(diff(t))); subplot(2,1,2); loglog(f,p,fn,pn,fnp,pnp,fnn,pnn) set(gca, 'ylim' ,yl); xlabel( 'f ' ); ylabel( '\\phi_x ' ); axis tight print -dpng -r250 doc/AliasingFixed Finite record length The last foible is finite record length. The difference between frequency bins is df = 1/T where T is the length of the record. Shorter records, the longer df . To see the effect make a second peak very close to the first, but far enough apart that they can be distinguished in the full time-series. om = th.f*2*pi; om_peak = find(om0.4*2*pi); f0=om(om_peak(1))/2/pi; om_peak=om(om_peak(1)+3); % chose a peak 4 bins higher. ff = om_peak/2/pi; 1./abs(ff-f0) ph = rand(1,1)*2*pi; sig.x_peak2 = A_peak.*cos(om_peak*sig.t+ph)/10; x = sig.x_peak2+sig.x_peak+sig.x_w; jmkfigure(3,2,0.7);clf subplot(3,1,1); plot(x); hold on ; plot(x(1:1500), 'b' ); plot(x(1:1000), 'r' ); ylabel( 'x ' ); subplot(3,1, ); sty = { 'k-' , 'b-' , 'r-' }; N = length(x); for i=1:3; =rawPowerSpec(x(1:N-0.35*N*(i-1)),dt); loglog(f,p,sty{i}); hold on ; h(i)=loglog(f,p, 'o' , 'markersi' ,6, 'col' ,sty{i}(1)); set(gca, 'xlim' , ); df(i)=median(diff(f)); nn(i) = length(x(1:N-0.25*N*(i-1))); str{i}=sprintf( 'N=%d; \\\\Delta f=%1.2e ' ,nn(i),df(i)); end ; xlabel( 'f ' ); ylabel( '\\phi_x ' ); legend(h(1:3),str{1:3}) title(sprintf( 'Peaks separated by \\\\Delta f=%1.2e' ,abs(ff-f0))); print -dpng -r400 doc/FiniteRecord.jpg ans = 66.6667 Here it is clear that fewer data points means poorer peak resolution. Once df rises above the peak separation they merge. A trick to get better frequency resolution is to zero-pad the data. However, in this case, even if we do this, we do not improve resolution between the peaks, just the centering of the peaks. jmkfigure(4,2,0.6);clf sty = { 'k-' , 'b-' , 'r-' }; N = length(x); for i=1:3; subplot3(1,3,i, 'xscale' , 'log' , 'yscale' , 'log' ); =rawPowerSpec(x(1:N-0.35*N*(i-1)),dt); loglog(f,p,sty{i}); xx = x(1:N-0.35*N*(i-1)); xxx = zeros(1,floor(5*N)); xxx(1:length(xx)) =xx; =rawPowerSpec(xxx,dt); h(i)=loglog(ff,pp,sty{i}); hold on ; h(i)=loglog(f,p,sty{i}); hold on ; % h(i)=loglog(f,p,'o','markersi',6,'col',sty{i}(1)); set(gca, 'xlim' , ); df(i)=median(diff(f)); nn(i) = length(x(1:N-0.25*N*(i-1))); str{i}=sprintf( 'N=%d' ,nn(i)); if i==1 ylabel( '\\phi_x ' ); end ; end ; %legend(h(1:3),str{1:3}) title( 'Zero-padded data' ); xlabel( 'f ' ); print -dpng -r400 doc/ZeroPad.jpg Note in this example how the poorly-resolved peaks (blue) have better peak localization than if we don't zero pad. However, this doesn't help separate peaks that aren't resolved in the unpadded timeseries. Spectral Leakage: Compute a raw fft of a strong peak with white noise. Note that we have put the steps from the previous tutorial into a function. dt = median(diff(sig.t)) =rawPowerSpec(sig.x_peak*10,dt); =rawPowerSpec(sig.x_w,dt); =rawPowerSpec(sig.x_peak*10+sig.x_w,dt); figure(1);clf semilogy(f,p,f,p_w,f,p_p); dt = 0.1000 Next Next we consider ways to reduce smearing of peaks in the spectra and ways to improve the signal-to-noise of spectral extimates. Published with MATLAB® 7.7 原文请访问：http://valdez.seos.uvic.ca/~jklymak/Phy580/html/FftFoibles.html
FFT first tutorial.. Contents Introduction Synthesized data Theoretical power spectra Introduction This tutorial deals with Fourier analysis and its use in interpreting data. The reader is assumed to be familiar with Fourier transforms, but not necessarily with their use in data characterization. Synthesized data Here we create a synthetic data set consisting of a red signal, an isolated peak, and white noise. First, the time axis is determined. Note this could be a spatial axis as well. This data has a sampling frequency of dt . We could use dt=1 , but this lends considerable confusion later when determining the amplitude of the spectrum. dt = 0.1; % seconds t = (0:2000)*0.1; We will have data with frequencies of 1/dt/2 to 1/1000; Note that we keep track of units as factors of 2*pi can easily be lost when dealing with the frequency domain. dom = 1/200; om = dom:dom:1/dt/2; % this is cycles/s=Hz. om = 2*pi*om; % this is rad/s. set the relative amplitudes: A_w = 0.3; A_r = 1; A_peak = 22; om_peak = find(om0.4*2*pi); % make sure it is not dead on a peak om_peak=om(om_peak(1))+0.3*median(diff(om)); White noise is a Gaussian process... x_w = A_w*randn(1,length(t)); figure(1); subplot(4,1,1); plot(t,x_w, 'col' , *0.5); A Red signal is the integral of a Gaussian process. %A = repmat(A_r*(1+0.9*randn(length(om),1)).*om'.^(-1),1,length(t)); %ph = repmat(rand(length(om),1)*2*pi,1,length(t)); %x_r = sum(A.*cos(om'*t+ph))*sqrt(dom); x_r = cumsum(randn(1,length(t)))*median(diff(t)); subplot(4,1,2); plot(t,x_r, 'r' ); A peak ph = rand(1,1)*2*pi; x_peak = A_peak.*cos(om_peak*t+ph)*sqrt(dom); subplot(4,1,3); plot(t,x_peak, 'b' ); Combine... x=x_peak+x_r+x_w; subplot(4,1,4); plot(t,x); set(datachildren(gcf), 'ylim' , *10); ylabel( 'x ' ); Theoretical power spectra The three wave forms plotted above have known power spectra. f = om/2/pi; df = median(diff(f)); % it is usual to plot spectra in terms of cycles/s. P_w = om*0+var(x_w)/max(f); P_r = f.^(-2); P_r = P_r./sum(P_r*df)*var(x_r); P_peak = om*0; in = find(om=om_peak);in =in(1) P_peak(in) = 1*var(x_peak)/df; % rename so we don't overwrite by accident. th.f=f; th.P_r=P_r; th.P_w=P_w; th.P_peak=P_peak; th.Px = th.P_r+th.P_w+th.P_peak; % sig.t =t; sig.x_w = x_w; sig.x_r = x_r; sig.x_peak = x_peak; sig.x = x; jmkfigure(2,2,0.4);clf loglog(th.f,th.P_w, 'col' , *0.5); hold on ; loglog(th.f,th.P_r, 'r' ); loglog(th.f,th.P_peak, 'bx' , 'markersi' ,10); loglog(th.f,th.Px, 'k' ); xlabel( 'f ' ); ylabel( 'P_x ' ); legend( 'x_w' , 'x_r' , 'x_{peak}' , 'Sum' ); title( 'Theoretical spectrum of synthetic signals' ); in = 82 The plot shown above is the ideal we are striving towards. With infinite samples of the processes the power spectra would look like what is given above. save SynthSig sig th if 0 warning off ; opt.outputDir= '/Users/jklymak/Sites/Phy580/html' ; publish( 'FftSynthesizeData' ,opt); end ; Published with MATLAB® 7.7 原文请访问： http://valdez.seos.uvic.ca/~jklymak/Phy580/html/FftSynthesizeData.html
标签: 谱分析
