i have results given by instrutor ,but i need MATLAB code for this
Assignment #1: estimation of correlation functions, PSD functions, and relationships for LTI systems inputs and outputs
Note: this document uses the same notation as in the course slides for Section 2, e.g., mx , mˆx , fyx( )l , fˆyx( )l , pyx( )l , Fyx(ej?), Pyx(ej?), and PyxW (ej?).
Note: equations in red correspond to additional information which does not appear in the slides.
Throughout the assignment, we use the following relationships describing two LTI systems:
, with impulse response g n( ) 5 (= dn- +1) 4 (dn-2)
k=-8
n
, with impulse response hn( )= ? ?? ?1 un( ).
k=-8 ? ?2
Each signal x n( ), d n( ), y n( ) is an observation (sequence, realization) of L samples from a wide sense stationary, mean ergodic and correlation ergodic random process X .
The signal x n( ) has an autocorrelation function of: f sdxx( )l = x2 ( )l + mx 2 with sx2 =1 and mx = 1 4 (although sometimes we will usemx = 0).
Question 1
a) Find the theoretical expressions for fyx()l ,fdx()l ,fyy()l .
b) Use the following code to generate an instance of the sequences x n( ), d n( ), y n( ) with L samples:
L=2^16; % nb. of samples for each realization
g=[0,5,4]; % FIR filter coeffs hb=[1]; % IIR filter numerator ha=[1,-1/2]; % IIR filter denominator
extra_samples_to_remove_transients = max(length(g),5/(1-abs(ha(2))))-1;
sigma_x_sq=1; % variance of input x mx=1/4; % mean of input x
x=sqrt(sigma_x_sq)*randn(1,L+extra_samples_to_remove_transients)+mx; d=filter(g,1,x); y=filter(hb,ha,x);
% remove transients
x=x(extra_samples_to_remove_transients+1:end); d=d(extra_samples_to_remove_transients+1:end); y=y(extra_samples_to_remove_transients+1:end);
Compute the unbiased and biased cross-correlation estimates fˆyx( )l and pyx( )l and plot them on the same figure as the theoretical fyx()l (e.g., three subplots). Produce two figures, using the two following scales for the “lag” axis: -10= =l 10 and - +L 1= =l L -1.
Note: for all the plots in the assignment, you must include meaningful titles, label for horizontal axis, and a grid.
Note: in this assignment you can use either the plot() or stem() functions, depending on which one produces the clearest plot for each case.
Note: If you use functions such as xcorr() instead of programming the details of the estimated correlation function, make sure that you use it with the correct parameters such that it gives the same results as our definitions (i.e., same scale, same normalization, etc.).
Repeat for fˆdx()l ,pdx()l ,fdx()l , and fˆyy( )l ,pyy()l ,fyy()l .
c) Generating N =100 sequences of y n( ), you can experimentally approximate the variance of the unbiased estimate fˆyy( )l with the following expression:
where fˆyyi ( )l represents the autocorrelation estimate measured from the i th sequence of y n( ).
Likewise, the variance of the biased autocorrelation estimate pyy()l can be experimentally measured, considering the known bias of this estimate:
?? ? l? ?2?2
var{pyy()l } = E ??(pyy()l - E p?? yy()l ??)2??? = E ??????pyy()l -???1- L ???fyy()l ??? ???? ˜ N1 ?iN=1 ????pyyi ()l -????1- Ll ????fyy()l ????
?
?
In principle, the range - +L 1= =l L -1 should be used for the plots, but you will need to discard the larger values of l near - +L 1 and L -1 (e.g. cut 1% on each side) because the measured variance at those lags becomes too large and prevents seeing the rest of the curves.
d) It is also possible to compare the measured var{fˆyy()l } and var{pyy()l } (previous question) with corresponding theoretical values:
varth {fˆyy()l } ˜ (L Ll )2 ??? n=- +?L-L1 1fyy( )n 2 + n=- + +L?- -L1 1ll f*yy(n - l) (fyy n + l)???? = (L -Ll )2fˆtmp()l
- ?
varth {pyy()l } ˜ L1 ??? n=- +?L-L1 1fyy( )n 2 + n=- + +L?- -L1 1l l f*yy(n - l) (fyy n + l)???? = L1fˆtmp()l - +L 1= =l L -1
?
For example a code like this can be included in your code to compute such functions: vartmp=ones(1,2*L-1)*sum(abs(phi_yy.*phi_yy)); for lag=0:L-1
vartmp(L+lag)=vartmp(L+lag)+sum(phi_yy(1:end-2*lag).*phi_yy(1+2*lag:end)); end
vartmp(1:L-1)=flip(vartmp(L+1:end));
var_phi_yy_hat_th=vartmp.*(L./(L-abs(-L+1:1:L-1)).^2); var_p_yy_th=vartmp/L;
However, these equations are only valid for mx = 0 (or they should be applied for autocovariance functions instead of autocorrelation functions), and for fairly small absolute lag l values. So here you must repeat the simulations of c) but using mx =0
. Plot var{fˆyy()l },var{pyy()l },varth {fˆyy()l },varth {pyy()l } on the same figure (four subplots) for the following ranges of lags (i.e. two figures):
• Using 15% of lag values (i.e., smaller absolute values around l =0); • From -200= =l 200.
Question 2
Note: for this question on PSD estimation, you must not use Matlab/Octave functions such as pwelch(), periodogram(), or cpsd() to compute the PSD quantities. To better understand the computations, use the basic fft() function and the equations for the PSD estimates. Note: frequency domain signals in this question are written with the discrete time Fourier transform (DTFT) notation (e.g. Fxx(ej?) -p ? p= ), but it is understood that in practice they are discretized/computed with fast Fourier transforms (FFTs):
Fxx( )k 0 = =k Nfft -1, ? pk = k2Nfft , where Nfft is the size of the FFT (after zeropadding if zero-padding is used). FFT results like Fxx( )k can be plotted to approximate DTFTs, with proper scaling of frequency axis and “wrap-around” or “fftshift” processing so that frequencies p ? p= 2 are shown as -p ?= 0.
a) Find the theoretical expressions for Fyx(ej?), Fyy(ej?), for the cases with mx =1 4 and mx =0.
b) For the cases with mx = 1 4 and mx = 0, generate x n( ) and y n( ) sequences of L = 216 samples, compute the periodogram estimates Pyx(ej?)and Pyy(ej?), and plot Pyx(ej?) ,
Pyy(ej?), Fyx(ej?) , Fyy(ej?) in the same figure (four subplots). Note that for mx = 1 4 the functions Fyx(ej?), Fyy(ej?) have infinite Dirac components which cannot be plotted
(just ignore them for the plots).
You will see that for the case mx = 1 4 the DC component (frequency zero) in Pyx(ej?) and
Pyy(ej?) completely hides the rest of the plot. But if you zoom the results with mx = 1 4 so that small amplitudes are displayed (use another figure for that), you will then see that except for the DC component the results are the same as for mx =0.
Therefore, it is often convenient in practice to remove the DC component of a signal before computing its Fourier transform or PSD estimate. For the rest of this assignment, generate signals with mx = 0.
c) From the previous part b), you should notice that the periodogram estimates Pyx(ej?), Pyy(ej?) have a lot of fluctuations: this is the effect of the large variance of these estimators, which does not decay with the length L of the observed signals. In this part, you are asked to verify this fact.
Generating N =100 sequences of y n( ), you can experimentally measure the variance of the periodogram estimate Pyy(ej?) with the following expression: var ( ) ( ) ( ) 2 ,
=1
where mˆP = N1 ?iN=1 P eyyi ( j?), and Pyyi (ej?) represents the periodogram estimate measured from the i th sequence of y n( ). Using different values of L (ex. L = 212 and L = 216 ), verify that the variance of the periodogram estimate Pyy(ej?) does not significantly decay as the number of observed samples L is increased.
Also compare (by plotting on the same figure) with the theoretical variance of the periodogram estimate:
Var P{ yy(ej?)} = E P[( yy(ej?) - E P[ yy(ej?)]) ]2 ˜Fyy2 (ej?).
d) With L = 216, M = 128, K = +1 ??(L M M- )/( / 2)?? = 1023 (case with 50% M /2 samples overlap between windows), and using a (normalized) Hamming window, find the Welch estimate PyyW (ej?). Plot PyyW (ej?) and the theoretical Fyy(ej?) on the same figure (two subplots).
e) Generating N =100 sequences of y n( ) with the same parameters L = 216, M = 128, K = +1 ??(L - M)/(M / 2)?? = 1023 , you can experimentally measure the variance of the Welch estimate PyyW (ej?) with the following expression: var ( ) ( ) ( ) 2 ,
=1
where mˆP = N1 ?iN=1 PyyW i, (ej?), and PyyW i, (ej?) represents the Welch estimate measured from the i th sequence of y n( ).
Also compare (by plotting on the same figure) with the theoretical variance of the Welch estimate for 50% overlap:
Var P e{ yyW ( j?)} ˜ 89K Fyy2 (ej?).
Note the difference between the variance of the Welch estimate var{P eyyW ( j?)} and the previous variance of the periodogram estimate var{P eyy( j?)}.
f) With L = 216, M = 128, K = +1 2(L M M- )/ = 1023 (case with 50% overlap between windows), and using a (normalized) Hamming window, compute an estimate PyxW (ej?) of the cross-PSD between x n( ) and y n( ), compute estimates of the auto-PSDs PxxW (ej?),
PyyW (ej?), and use them to compute an estimated coherence function
?Wyy(ej?)= PyxW (ej?) PxxW (ej?)PyyW (ej?). Compare ?Wyy(ej?) with the theoretical coherence magnitude ?yx(ej?) by plotting them on the same figure.