-- Blackman-Tukey FFT MEM Blackman-Tukey MEM MEM MEM MEM Singular Spectrum Analysis Multi-Taper Method (Matlab pmtm) 3... y(t) (Fourier transform) t=

Size: px

Start display at page:

Download "-- Blackman-Tukey FFT MEM Blackman-Tukey MEM MEM MEM MEM Singular Spectrum Analysis Multi-Taper Method (Matlab pmtm) 3... y(t) (Fourier transform) t="

さいぞうすみい
5 years ago
Views:

1 Matlab MEM MEM MEM Matlab MEM y(t) (Blackman-Tukey) y(t) (periodgram) (Maximum Entropy Method: MEM)

2 -- Blackman-Tukey FFT MEM Blackman-Tukey MEM MEM MEM MEM Singular Spectrum Analysis Multi-Taper Method (Matlab pmtm) 3... y(t) (Fourier transform) t= t= ( π ) Yk ( ) = yt ( )exp i kt/ = yt ( )exp( i π ft t); f = k/( t), k= 0,..., k k Y / t, / t,... k k f = /( t) (yquist frequency) 0.5 (3.) exp( i π( k) / ) = cos( π( k/ )) isin( π( k/ )) = cos( π) cos( πk/ ) i cos( π) sin( πk/ ) = cos( πk/ ) + isin( πk/ ) = exp( + i π k/ ) /( t) < fk < / t f < fk < f f = /( t) (fundamental frequency) (3.) (inverse Fourier transform) (3.)

3 -3- iπ kn yt ( ) = Y( k)exp, t=,...,, f= k = 0 t (/) Matlab fft (time domain) (frequency domain) (periodogram) (Power Spectrum Density, PSD), S yy Syy ( fk ) = Y( k), (3.) fk, Yk (3.) f < fk < f Syy ( fk ) = Y( k), f < fk f (3.3) (two-sided Power Spectrum Density) Y( f) = Y( f ) * S ( f ) = S ( f ) yy k yy k (one-sided Power Spectrum Density),...,, Gyy ( fk ) = Yk, k = 0, Gyy ( 0 ) = Y0,, Gyy ( f /) = Gyy ( f /) = Y /, k / (3.4) / yn = Gyy (3.5) n= n= 0 S yy y S xy

4 -4- Syy ( fk ) = Y( fk ), 0 fk f (3.6) yy ( k ) ( k ), 0 k S f = Y f f f (3.7) (Parseval s energy conservation theorem) yn = Syy (3.8) n= n= % 50% yn, n = 0, i yt () = Yk ( )exp k = 0 π kt

5 -5- n= 0 yt () * i kn iπkn yt () = Yk ( )exp Yk ( )exp π t= 0 n= 0 k= k= * ( n ) f = /( t) yt ( ) = Yk ( ) = Yk ( ) n= 0 t= 0 k= 0 k= 0 (3.8) yt () = Y( k) = S = = t= 0 k= 0 k= 0 k (Wiener-Khinchin relation) 3. ( ) (/ ) ρτ = x ( t ) x ( t+ τ ) t= 0 ( )e = e xtxt ( ) ( + τ) ( i πkτ / ) ( i πkτ / ) ρτ τ= 0 τ= 0 t= 0 = xt xt+ τe e t= 0 τ = 0 i π ( k) t/ = xe t { X( k) } t= 0 = X( k) X( k) = X( k) i πk( t+ t)/ i πk( t)/

6 Welch (Welch / m m/ m (smoothing window) (data window) (time window) (taper) (spectral window)

7 (spectral leakage) (taper) ht () yt () Y '( k) = t h( t) y( t) exp( i π kt f), k =,...,, f = /( t) t= Y( k) yt ( ) = f Yl ( ) exp( i π lt f) t= l= l= Y '( k) = t h( t) f Y( l) exp( i πlt f ) exp( i πkt f ) Y'( k) = f Y( l) t h( t)exp( i π ( k l) t f) l= t= = f Y() l H( k l) l= H h(t) t= Hk ( ) = t ht ( ) exp( i π kt f), (3.9) S'( k) = Y '( k) (3.9) S'( k) (99X, p. 86) Percival and Walden (993, p. 07) { } (3.9) f (3.0) f S' = H( f f ') S( f ') df ' (Hamming taper)

8 -8- (Hanning taper) w ( n t) = π =,... [.0 cos( ( n 0.5) / )], n, ( ( n 0.5) / ), n, w ( n t) = cos π =,... Hamming taper Hanning taper 0% cos cos0 cos0 hamming, hanning taper 3.3 Hanning( )Hamming ( )cos0( ) yt ( ) = cos( πt/ ) + 0.cos( πt/ 9) 00 boxcar tapercos0, hamming 95% X boxcar cos0 hamming taper hamming taper 3.4 yt ( ) = cos( πt/ 50) + cos( πt/ p), p= 3 + (7 3)( t ) / 99 95% main lobe

9 hamming taper boxcar cos0( )hamming ( ) 0 log0 8 95%

10 FFT direct method rectangle taper (977) 0% cosine Hamming window (989)Percival and Walden (993) K (zero-padding) 3 5

11 -- S f = Y( f ) = y exp( i π f n t), (3.) + K yy ( k ) k n k K n= 0 ( + ) f = k/(( + K) t), f f f k k (direct spectrum estimator) (modified periodogram) (directly smoothed direct spectrum estimator) (Percival and Walden, 993) 00 00, 50, 33, 5, % side lobe (ripple) hamming

12 % Matlab Matlab psdpwelch tsrs x=(::00) ; tsrs=sin(x**pi/0); subplot(,,); plot(x,tsrs) % tlng=length(tsrs); % nfft=56; %fft nfft-tlng avfnc=ones(5,); % avfnc=avfnc./sum(avfnc); % taper=hamming(tlng); % hamming taper tsrs_tapered=tsrs.*taper; %

13 -3- tsrs_tapered(tlng+:nfft)=0.; % famp=fft(tsrs_tapered); % famp=famp(:nfft/+)/tlng; % ds_psd=abs(famp.^); % dsm_psd=conv(ds_psd,avfnc); % nav=(length(avfnc)-)/; dsm_psd = dsm_psd (nav+:end-nav); % frq=(0:nfft/)./nfft; % subplot(,,); plot(frq,dsm_psd, frq,ds_psd) n n R I χ = ( Y ) ( ) ˆ i + Yi Y = ns S n S i= Ŝ P χ n α α < χ < χ n = α χ = ns ˆ / S P nsˆ nsˆ < S < = α α α Fχ,n Fχ,n (3.) F χ ICDF α α Fx,n Fx,n, Fx ( 0.5,n) Fx ( 0.5,n) (3.3)

14 -4- (equivalent degree of freedom, EDOF) Q s ( f ) = s( f ) Q( f f ') df ' Q (smoothing window) Q( f ) df = (equivalent band width) σ n n i= = (/ ) = / n Q ( f ) df, B e, B e =, (3.4) Q ( f) df e = Be (3.)(3.3) Tt () h spectral window (994) Percival and Walden (994) p. 38 spectral window von Stroch and Zwiers (999)

15 -5- ( ) e = + Be h K e = + ( Be ) h + K (3.5) (3.5) y' ( t) = y( t) ry( t ) y' pre-whitening psd_ex Matlab function [ds_psd,frq]= psd_ex(tsrs,nfft,ctaper,avfnc) % function [ds_psd,frq]= psd_ex(tsrs,nfft,ctaper,avfnc) % Power Spectrum Density Matlab Script % % sz> % % tsrs: (sz, sz) % nfft: fft % ctaper: taper ('hamming','hanning' ) % avfnc:

16 -6- % % ds_psd: nfft/+ % frq: nfft/+ % todo % PSD [sz,sz]=size(tsrs); if (sz==); tsrs=tsrs(:); end %tsrs sz* if (sz>nfft) error('*psd_ex.m* sz>nfft') %nfft tlng end % avfnc=avfnc./sum(avfnc); % str=sprintf('%s(%d)',ctaper,sz); %taper hamming(00) taper=eval(str); % taper for n=:sz tsrs(:,n)=tsrs(:,n).*taper; % end tsrs(sz+:nfft,:)=0.; % famp=fft(tsrs); % famp=famp(:nfft/+,:)/sz; % d_psd=abs(famp.^); % (d: direct) for n=:sz ds_psd(:,n)=conv(d_psd(:,n),avfnc); % (ds: direct-smoothed) end nav=(length(avfnc)-)/; ds_psd=ds_psd (nav+:end-nav,:); % frq=(0:nfft/)./nfft; % if (nargout==0) plot(frq,ds_psd); xlabel('frequency'); ylabel('psd') end return for npnl=:4 tlng=00; nsmpl=000; if (npnl==) nfft=00; taper=boxcar(tlng); ctaper='boxcar'; avfnc=[ ]; cttl='no padding, boxcar, 5-running'

17 -7- elseif (npnl==) nfft=00; taper=hamming(tlng); ctaper='hamming'; avfnc=[ ]; cttl='no padding, hamming, 5-running' elseif (npnl==3) nfft=00; taper=boxcar(tlng); ctaper='boxcar'; avfnc=[ ]; cttl='double padding, boxcar, 5-running' elseif (npnl==4) nfft=00; taper=hamming(tlng); ctaper='hamming'; avfnc=conv([ ],[ ]); cttl='double padding, hamming, double --' end tsrs=randn(tlng,nsmpl); [ds_psd,frq]=psd_ex(tsrs,nfft,ctaper,avfnc); % Monte-Carlo 90% srt_psd=sort(ds_psd,); clear c; c(:,)=log0(srt_psd(:,nsmpl*0.05)); c(:,)=log0(srt_psd(:,nsmpl*0.95)); avfnc=avfnc./sum(avfnc); % frequency domain c0=srt_psd(:,nsmpl*0.5); edof=(+(/sum(avfnc.^)-)*tlng/nfft*mean(taper)) clear ce; ce(:,)=log0(c0*chiinv_nint(0.050,*edof)/chiinv_nint(0.50,*edof)); ce(:,)=log0(c0*chiinv_nint(0.950,*edof)/chiinv_nint(0.50,*edof)); subplot(,,npnl) plot(frq,c,'b',frq,ce,'r'); title(cttl); xlabel('frq'); ylabel('psd'); if (npnl==); legend('monte-carlo','','theory'); end end

18 % 3.4. MEM MEM MEM MEM Press (993) MEM (peak splitting) (false peak) MEM Percival and Walden (99) MEM AR

19 MEM MEM yt ( ) + φ() yt ( ) φ( pyt ) ( p) + σwε( t) = 0 { yt ( ) + φ() yt ( ) φ( pyt ) ( p) } exp( i πkτ / ) = σw ε( t)exp( i πkτ / ) (3.6) t= 0 t= 0 ( ) ( ) ( ) y( t τ)exp i πkt/ = y( t τ)exp i πk( t τ)/ exp i πkτ / t= 0 t= 0 ( π τ ) = Yk ( ) exp i k / (3.6) { φ ( π ) φ ( π ) φ ( π )} = Yk ( ) + ()exp i k/ + ()exp i4 k/ ( p)exp i p k/ p = Yk ( ) + φ()exp( i πkτ / ) τ = p + φ π τ = σw τ = Yk ( ) ()exp( i k / ) MEM Yk ( ) σw = p + φ() exp( i πkτ / ) τ = (3.7) MEM (all pole) MA ARMA Matlab MEM Matlab MEM pmem, pyulear, pburg pyulear pburg Yule-Walker Burg pmem tlng=00;

20 -0- rtsrs=randn(tlng,); fs=; order=0; nfft=00; %fs [pm,frq]=pmem(rtsrs,order,nfft,fs); [pm,frq]=pburg(rtsrs,order,nfft,fs); [pm3,frq]=pyulear(rtsrs,order,nfft,fs); plot(frq,pm,frq,pm,frq,pm3) MEM MEM MEM MEM %----- <phase randamization > tlng=00; nfft=00; nh=nfft/; nhp=nfft/+; nsmpl=000; frqx=([::nhp]'-)/nhp*0.5; frqc=0.; fwid=0.0; amp=.0; spct_gvn=amp*exp(-(frqx-frqc).^./fwid^); frqc=0.05; fwid=0.0; amp=0.7; spct_gvn=amp*exp(-(frqx-frqc).^./fwid^)+spct_gvn; frqc3=0.5; fwid3=0.0; amp3=0.0; spct_gvn=amp3*exp(-(frqx-frqc3).^./fwid3^)+spct_gvn; nlvl=0.0; spct_gvn=spct_gvn+nlvl; subplot(3,,); plot(frqx,log0(spct_gvn)); title('given spectrum') xlabel('frq'); ylabel('log0(psd)') famp=sqrt(spct_gvn); famp=famp(:)*ones(,nsmpl); fphs=rand(length(frqx),nsmpl)**pi; fcmp=famp.*exp(i*fphs); fftamp=zeros(nfft,nsmpl); fftamp(:nhp,:)=fcmp; fftamp(nhp+:end,:)=conj(flipud(fcmp(:nh,:))); rtsrs=real(ifft(fftamp)); rtsrs=rtsrs(:tlng,:); x=[::tlng]'; subplot(3,,); plot(x,rtsrs(:,)); title('sample time series') xlabel('time'); ylabel('y') %----- </phase randamization > %----- < > ctaper='boxcar'; nfft=00; avfnc=[ ]; [ds_psd,frq_ds]=psd_ex(rtsrs,nfft,ctaper,avfnc); ctaper='hamming'; nfft=00; avfnc=[ ]; [ds_psd,frq_ds]=psd_ex(rtsrs,nfft,ctaper,avfnc); %----- </ > %----- <MEM > fs=.0; nfft=56; clear pmy pmb; order=33; for n=:nsmpl [pmb(:,n),frq]=pyulear(rtsrs(:,n),order,nfft,fs); end

21 -- order=0; for n=:nsmpl [pmb(:,n),frq]=pyulear(rtsrs(:,n),order,nfft,fs); end %----- </MEM > clear pmyest pmbest pdsest pdsest; srt_dsp=sort(ds_psd,)*tlng; pdsest(:,)=srt_dsp(:,0.05*nsmpl); pdsest(:,)=srt_dsp(:,0.5*nsmpl); pdsest(:,3)=srt_dsp(:,0.95*nsmpl); subplot(3,,3); plot(frq_ds,log0(pdsest)); title('direct-psd, double-zero, --, boxcar') xlabel('frq'); ylabel('log0(psd)') srt_dsp=sort(ds_psd,)*tlng; pdsest(:,)=srt_dsp(:,0.05*nsmpl); pdsest(:,)=srt_dsp(:,0.5*nsmpl); pdsest(:,3)=srt_dsp(:,0.95*nsmpl); subplot(3,,4); plot(frq_ds,log0(pdsest)) title('direct-psd, no-zero, --, hamming') xlabel('frq'); ylabel('log0(psd)') srt_pmb=sort(pmb,); pmbest(:,)=srt_pmb(:,0.05*nsmpl); pmbest(:,)=srt_pmb(:,0.5*nsmpl); pmbest(:,3)=srt_pmb(:,0.95*nsmpl); subplot(3,,5); plot(frq,log0(pmbest)) title('pburg, 33-order'); xlabel('frq'); ylabel('log0(psd)') srt_pmb=sort(pmb,); pmbest(:,)=srt_pmb(:,0.05*nsmpl); pmbest(:,)=srt_pmb(:,0.5*nsmpl); pmbest(:,3)=srt_pmb(:,0.95*nsmpl); subplot(3,,6); plot(frq,log0(pmbest)) title('pburg, 0-order'); xlabel('frq'); ylabel('log0(psd)') legend('5%','50%','95%')

22 MEM

23 -3-5% 5% =0.9 9%.5 5% 5% 994:, p. 99. Percival D. B. and A. T. Walden, 993: Spectral analysis for physical applications: Multaper and conventioinal univariate techniques. Cambridge University Press, pp Press, W. H., B. P. Flannery, S. A. Teukolsky, W. T. Vetterling, 993: umerical Recipes in C [ ],, pp von Stroch, H. and F. W. Zwiers, 999: Statistical analysis in climate research. Cambridge University Press, pp. 483 ISB:

[1] 1.1 x(t) t x(t + n ) = x(t) (n = 1,, 3, ) { x(t) : : 1 [ /, /] 1 x(t) = a + a 1 cos πt + a cos 4πt + + a n cos nπt + + b 1 sin πt + b sin 4πt = a

13/7/1 II ( / A: ) (1) 1 [] (, ) ( ) ( ) ( ) etc. etc. 1. 1 [1] 1.1 x(t) t x(t + n ) = x(t) (n = 1,, 3, ) { x(t) : : 1 [ /, /] 1 x(t) = a + a 1 cos πt + a cos 4πt + + a n cos nπt + + b 1 sin πt + b sin