Spectral Estimation Overview Introduction Periodogram R x (e j - PowerPoint PPT Presentation

Spectral Estimation Overview Introduction ∞ • Periodogram � R x (e jω ) � r x ( ℓ )e − jωℓ • Bias, variance, and distribution l = −∞ • Blackman-Tukey Method • Most stationary random processes have continuous spectra • Welch-Bartlett Method • If we have time, will discuss line spectra at end of term • Others • Recall the definition of PSD above • The estimation problem is to find a ˆ R x (e jω ) given only a finite data record { x ( n ) } N − 1 0 • If the autocorrelation can be estimated with great accuracy at long lags, can treat as a deterministic problem • With short records (or equivalently, local stationarity), is more difficult J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 1 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 2 Example 1: Data Windowed Autocorrelation Periodogram Solve for the expected value of the biased estimate of autocorrelation 2 � N − 1 � R x (e jω ) � 1 = 1 � � � 2 applied to a windowed data segment. How should w ( n ) be scaled such ˆ � v ( n )e − jωn � � V (e jω ) � � � N N � � that the estimate is unbiased at ℓ = 0 . � n =0 � ∞ where v ( n ) � x ( n ) w ( n ) . r x ( ℓ ) = 1 � x ( n + | ℓ | ) w ( n + | ℓ | ) x ∗ ( n ) w ∗ ( n ) ˆ N • If w ( n ) = c (a constant), is called the Periodogram n = N −| ℓ |− 1 • If w ( n ) is not constant, is called the Modified Periodogram and w ( n ) is called the data window • Can be estimated quickly using the FFT • Is related to the biased autocorrelation estimate we discussed last time N − 1 ˆ � r v ( ℓ )e − jωℓ R x (e jω ) = ˆ ℓ = − ( N − 1) J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 3 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 4

Example 2: Periodogram and Biased Autocorrelation Example 2: Workspace Let x ( n ) be a finite duration sequence, 0 to N − 1 . Show that the biased estimate of autocorrelation is given by r x ( ℓ ) = 1 ˆ N x ( ℓ ) ∗ x ( − ℓ ) and that the DTFT of ˆ r x ( ℓ ) is equal to the Periodogram of x ( n ) when w ( n ) = 1 . J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 5 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 6 Periodogram Example 3: Periodogram as a Filter Bank 2 2 � N − 1 � � N − 1 � R x (e jω ) � 1 = 1 � 2 = F { ˆ R x (e jω ) � 1 � � � � ˆ � ˆ � v ( n )e − jωn � V (e jω ) � � v ( n )e − jωn r v ( ℓ ) } � � � � N N N � � � � � n =0 � � n =0 � • Remember that we must use the biased estimate of • How is the periodogram a filter bank? autocorrelation – Otherwise the estimated PSD may be negative at some • What is the transfer function of the filter? frequencies • How good are the filters? How close are they to ideal filters? � 2 = F { ˆ 1 � � � V (e jω ) r v ( ℓ ) } no longer holds – The equality N • This relationship to ˆ r ( ℓ ) means we can use the FFT to calculate ˆ r ( ℓ ) very efficiently – Take FFT of signal, magnitude square, and inverse FFT – Requires O ( N log N ) operations instead of O ( N 2 ) ! – Must take care to zero pad sufficiently J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 7 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 8

Example 3: Work Space Example 4: Peridogram Generate the periodogram for a signal from a known LTI system with two sets of complex conjugate poles close to one another and near the unit circle. Repeat for various signal lengths. J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 9 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 10 Example 4: MATLAB Code Example 4: Pole Zero Map of ARMA Process M = 500; % Padding to eliminate transient N = [25,50,100,250,500,1000]; % Length of signal segment 1 NA = 500; % No. averages cb = 95; % Confidence Bands NZ = 1024; 0.5 b = poly([-0.8,0.97*exp(j *pi/4),0.97*exp(-j *pi/4),... 0.97*exp(j *pi/6),0.97*exp(-j *pi/6)]); % Numerator a = poly([ 0.8,0.95*exp(j*3*pi/4),0.95*exp(-j*3*pi/4),... 0.95*exp(j*2.5*pi/4),0.95*exp(-j*2.5*pi/4)]); % Denominator 0 b = b*sum(a)/sum(b); % Scale DC gain to 1 −0.5 −1 −1 −0.5 0 0.5 1 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 11 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 12

Example 4: MATLAB Code Example 4: Signal figure Length:100 h = circle; 300 z = roots(b); p = roots(a); hold on; h2 = plot(real(z),imag(z),’bo’,real(p),imag(p),’rx’); 200 hold off; axis square; xlim([-1.1 1.1]); 100 ylim([-1.1 1.1]); axislines; Signal box off; 0 −100 −200 0 10 20 30 40 50 60 70 80 90 100 Sample Index J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 13 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 14 Example 4: MATLAB Code Example 4: True PSD n = 100; 5 x 10 w = randn(M+n,1); True PSD x = filter(b,a,w); % System with known PSD 3 nx = length(x); x = x(nx-n+1:nx); % Eliminate start-up transient (make stationary) 2 PSD figure; k = 0:n-1; 1 h = stem(k,x); set(h,’Marker’,’none’); 0 set(h,’LineWidth’,0.5); 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 box off; ylabel(’Signal’); xlabel(’Sample Index’); title(sprintf(’Length:%d’,n)); xlim([0 n]); ylim([min(x) max(x)]); PSD 0 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Frequency (cycles/sample) J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 15 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 16

Example 4: MATLAB Code Example 4: Periodogram Estimate [R,w] = freqz(b,a,NZ); 5 x 10 R2 = abs(R).^2; N:25 NA:1000 Confidence Bands:95% f = w/(2*pi); 4 Single Estimate subplot(2,1,1); h = plot(f,R2,’r’); True PSD 3 set(h,’LineWidth’,1.2); PSD Confidence Bands box off; 2 Average Estimate ylabel(’PSD’); title(’True PSD’); 1 xlim([0 0.5]); 0 ylim([0 max(R2)*1.2]); 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 subplot(2,1,2); h = semilogy(f,R2,’r’); set(h,’LineWidth’,1.2); box off; ylabel(’PSD’); xlim([0 0.5]); PSD ylim([0.2*min(R2) max(R2)*5]); 0 xlabel(’Frequency (cycles/sample)’); 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Frequency (cycles/sample) J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 17 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 18 Example 4: Periodogram Estimate Example 4: Periodogram Estimate 5 5 x 10 x 10 N:50 NA:1000 Confidence Bands:95% N:100 NA:1000 Confidence Bands:95% 4 4 Single Estimate Single Estimate True PSD True PSD 3 3 Confidence Bands Confidence Bands PSD PSD 2 2 Average Estimate Average Estimate 1 1 0 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 PSD PSD 0 0 10 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Frequency (cycles/sample) Frequency (cycles/sample) J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 19 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 20

Example 4: Periodogram Estimate Example 4: Periodogram Estimate 5 5 x 10 x 10 N:250 NA:1000 Confidence Bands:95% N:500 NA:1000 Confidence Bands:95% 4 4 Single Estimate Single Estimate True PSD True PSD 3 3 PSD Confidence Bands PSD Confidence Bands 2 2 Average Estimate Average Estimate 1 1 0 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 PSD PSD 0 0 10 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Frequency (cycles/sample) Frequency (cycles/sample) J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 21 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 22 Example 4: Periodogram Estimate Example 4: MATLAB Code for c1 = 1:length(N), 5 x 10 nx = M+N(c1); N:1000 NA:1000 Confidence Bands:95% Rh = zeros(NA,NZ/2+1); 4 Single Estimate kh = 1:NZ/2+1; True PSD fh = (kh-1)/NZ; 3 for c2 = 1:NA, Confidence Bands PSD 2 w = randn(M+N(c1),1); Average Estimate x = filter(b,a,w); % System with known PSD 1 x = x(nx-N(c1)+1:nx); % Eliminate start-up transient (make stationary) rh = (1/N(c1))*abs(fft(x,NZ)).^2; 0 Rh(c2,:) = rh(kh).’; 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 end; Rha = mean(Rh); % Average Rhu = prctile(Rh,100-(100-cb)/2); % Upper confidence band Rhl = prctile(Rh, (100-cb)/2); % Lower confidence band PSD 0 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Frequency (cycles/sample) J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 23 J. McNames Portland State University ECE 538/638 Spectral Estimation Ver. 1.14 24