Joint Signal Analysis Overview Introduction • Mostly we have focused on estimating statistical properties of a • Cross-correlation single univariate signal • Cross Power Spectrum – Autocorrelation function (ACF) • Examples – Partial autocorrelation function (PACF) – Power spectral density • In many applications we have two or more signals, x ( n ) and y ( n ) • Would like to say something about how they are related J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 1 J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 2 System Identification Joint Signal Analysis • Joint signal analysis characterizes the relationship between a pair w ( n ) G ( z ) of signals – We will focus on nonparametric methods – LTI systems x ( n ) y ( n ) H ( z ) Σ – Only 2 signals • Joint signal analysis is related to system identification • We have already discussed many of the possible properties • The goal of system identification is to build a model – Normalized cross-correlation, aka cross-correlation function (CCF)) • That is, estimate H ( z ) and G ( z ) , given x ( n ) and y ( n ) – Cross-power spectral density (CPSD) – Parametric, though order may be estimated – Coherence – Mostly LTI systems – Some methods for MIMO systems J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 3 J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 4
Joint Signal Analysis Causality • Strictly speaking, given only realizations of x ( n ) and y ( n ) we w ( n ) G ( z ) cannot determine whether one signal caused the other or not • This is a fundamental idea rooted in statistics x ( n ) y ( n ) H ( z ) Σ • However, under certain assumptions causality can be determined • For example, if we assume that any system that relates x ( n ) to • If we assume that x ( n ) and y ( n ) are related by an LTI system y ( n ) or vice versa is causal, then we may be able to get a sense of H ( z ) , then we can estimate the properties of H ( z ) given which caused which from the estimated cross-correlation realizations of x ( n ) and y ( n ) – Magnitude – Phase – Group delay – Impulse response J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 5 J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 6 Normalized Cross-Correlation Defined Estimated Cross-Covariance Normalized Cross-Correlation , also known as the Cross-correlation N − 1 −| ℓ | 1 Function (CCF) is defined for a WSS signal as � � ˆ γ u ( ℓ ) [ y ( n + | ℓ | ) − ˆ µ y ] [ x ( n ) − ˆ µ x ] N − | ℓ | γ yx ( ℓ ) = γ yx ( ℓ ) n =0 ρ yx ( ℓ ) = N − | ℓ | � σ y σ x γ y (0) γ x (0) � γ b ( ℓ ) ˆ ˆ γ u ( ℓ ) N where γ yx ( ℓ ) is the cross-covariance of y ( n ) and x ( n ) , N − 1 −| ℓ | 1 � = [ y ( n + | ℓ | ) − ˆ µ y ] [ x ( n ) − ˆ µ x ] E [( y ( n + ℓ ) − µ y )( x ( n ) − µ x ) ∗ ] γ yx ( ℓ ) = N n =0 for | ℓ | < N . ˆ γ u ( ℓ ) = ˆ γ b ( ℓ ) = 0 for | ℓ | ≥ N . • There are similar trade-offs between biased and unbiased estimates of the CCF as there were for the ACF • As with ACF, we generally assume the means are zero to simplify calculations and eliminate DC artifacts • The bias caused by estimating the means is usually negligible J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 7 J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 8
Estimated Cross-Covariance: Biased versus Unbiased Estimated Cross-Correlation Functions If µ y = µ x = 0 , The natural estimators of the CCF are ρ b ( ℓ ) � ˆ γ b ( ℓ ) ρ u ( ℓ ) � ˆ γ u ( ℓ ) N − 1 −| ℓ | ˆ ˆ 1 � σ x ˆ ˆ σ y σ x ˆ ˆ σ y γ u ( ℓ ) ˆ = y ( n + | ℓ | ) x ( n ) N − | ℓ | n =0 • In general, it is not possible to obtain confidence intervals for the N − 1 −| ℓ | estimated CCF because the variance of the estimator depends on 1 � ˆ γ b ( ℓ ) = y ( n + | ℓ | ) x ( n ) the true CCF and signal ACF’s N n =0 • Instead, it is common practice to plot the confidence intervals of a • The true CCF is not positive definite and the cross-spectral density purely random process is not non-negative, in general, so this is no longer a concern • If N is large enough, the central limit theorem applies and • Nonetheless, the unbiased estimator contains excessive variance at var { ˆ ρ b ( ℓ ) } is approximately normal for all ℓ large lags and thus the biased estimate is generally preferred • In this case, we can use the Normal cdf to plot confidence intervals of an IID sequence � • These are proportional to ± var { ˆ ρ yx ( ℓ ) } • This is the same as it was for the ACF J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 9 J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 10 Example 1: CCF Estimation Example 2: MATLAB Code Estimate the CCF of the ARMA process investigated earlier for PSD L = 200; % Length of autocorrelation calculated N = 250; estimation. cl = 99; % Confidence level NP = norminv((1-cl/100)/2); % Find corresponding lower percentile of Normal 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 b = b*sum(a)/sum(b); % Set DC gain to 1 x = randn(N,1); y = filter(1,a,x); J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 11 J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 12
Example 2: Pole Zero Map of ARMA Process Example 2: MATLAB Code figure 1 h = circle; z = roots(b); p = roots(a); hold on; h2 = plot(real(z),imag(z),’bo’,real(p),imag(p),’rx’); 0.5 hold off; axis square; xlim([-1.1 1.1]); ylim([-1.1 1.1]); 0 box off; −0.5 −1 −1 −0.5 0 0.5 1 J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 13 J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 14 Example 2: Signals Example 2: MATLAB Code figure N=250 n = 0:min(200,N); subplot(2,1,1); 2 h = stem(n,x(n+1)); set(h(1),’MarkerFaceColor’,’b’); x(n) set(h(1),’MarkerSize’,2); 0 xlim([0 max(n)]); ylim([min(x) max(x)]); ylabel(’x(n)’); −2 title(sprintf(’N=%d’,N)); box off; 0 20 40 60 80 100 120 140 160 180 200 subplot(2,1,2); h = stem(n,y(n+1),’g’); set(h(1),’MarkerFaceColor’,’g’); set(h(1),’MarkerSize’,2); 5 xlim([0 max(n)]); ylim([min(y) max(y)]); y(n) 0 ylabel(’y(n)’); xlabel(’Sample (n)’); −5 box off; −10 0 20 40 60 80 100 120 140 160 180 200 Sample (n) J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 15 J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 16
Example 2: Biased CCF Example 2: MATLAB Code np = 2^nextpow2(2*N-1); % Figure out how much to pad the signal N=250 1 Y = fft(y,np); X = fft(x,np); 0.8 cp = ifft(Y.*conj(X)); cc = real([cp(np-(L-1:-1:0));cp(1:L+1)]); 0.6 cc = cc/(std(y)*std(x)); cc = cc/N; 0.4 l = (-L:L).’; 0.2 figure h = stem(l,cc); ρ (l) 0 set(h(1),’Marker’,’None’); xlim([-L L]); −0.2 ylim([-1 1]); hold on; −0.4 plot(xlim,NP*[1 1]/sqrt(N),’k:’,xlim,-NP*[1 1]/sqrt(N),’k:’); plot([ 0 L],[1 (N-L)/N],’g’,[ 0 L],-[1 (N-L)/N],’g’); −0.6 plot([-L 0],[(N-L)/N 1],’g’,[-L 0],-[(N-L)/N 1],’g’); hold off; −0.8 ylabel(’\rho(l)’); xlabel(’Lag (l)’); −1 title(sprintf(’N=%d’,N)); −200 −150 −100 −50 0 50 100 150 200 box off; Lag (l) J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 17 J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 18 Example 2: Unbiased CCF Example 2: MATLAB Code np = 2^nextpow2(2*N-1); % Figure out how much to pad the signal N=250 1 Y = fft(y,np); X = fft(x,np); 0.8 cp = ifft(Y.*conj(X)); cc = real([cp(np-(L-1:-1:0));cp(1:L+1)]); 0.6 cc = cc/(std(y)*std(x)); l = (-L:L).’; 0.4 cc = cc./(N-abs(l)); 0.2 figure h = stem(l,cc); ρ (l) 0 set(h(1),’Marker’,’None’); xlim([-L L]); −0.2 ylim([-1 1]); hold on; −0.4 plot(l,NP*1./sqrt(N-abs(l)),’k:’,l,-NP*1./sqrt(N-abs(l)),’k:’); hold off; −0.6 ylabel(’\rho(l)’); xlabel(’Lag (l)’); −0.8 title(sprintf(’N=%d’,N)); box off; −1 −200 −150 −100 −50 0 50 100 150 200 Lag (l) J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 19 J. McNames Portland State University ECE 538/638 Joint Signal Analysis Ver. 1.05 20
Recommend
More recommend
