Matlab Code for Computer Exercise C8.6 Monson H. Hayes Statistical DSP

In this post, we will discuss Matlab Code for Computer Exercise C8.6 Monson H. Hayes Statistical Digital Signal Processing book. The output of matlab code simulation is available at the end of this page.

Matlab Code for computer exercise C8.6 is shown below.

%% Part (a)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];
w=filter(b,1,randn(64,1)); %data generation with gaussian noise - 100 samples
[Pxx,F] = pwelch(w);

plot(F,10*log10(Pxx));hold on
title('power spectrum of w(n)');
xlabel('Frequency');
ylabel('psd in dB');
grid on;
    x=2.*cos(n.*w1+phi1)+2.*cos(n.*w2+phi2)+w';
[Pxx_x,F_x] = pwelch(x);
plot(F_x,10*log10(Pxx_x));
title('power spectrum of x(n)');
xlabel('Frequency');
ylabel('psd in dB');
legend('psd of w(n)','psd of x(n)');
grid on;
% end part (a)-------------------------------------------

%% Part (b)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];
w=filter(b,1,randn(64,1));
va = var(w);
disp('The variance is:'),disp(va);
[Pxx,F] = pwelch(w);
mag=abs(Pxx);
 x1=2.*cos(n.*w1+phi1);
 x2=2.*cos(n.*w2+phi2);
 [Pxx_x1,F_x1] = pwelch(x1);
 mag_x1=abs(Pxx_x1);
 [Pxx_x2,F_x2] = pwelch(x2);
 mag_x2=abs(Pxx_x2);
plot(F,10*log10(Pxx));hold on
plot(F_x1,10*log10(Pxx_x1));hold on
plot(F_x2,10*log10(Pxx_x2));hold on
title('power spectrum');
xlabel('Frequency');
ylabel('psd in dB');
legend('psd of w(n)','psd of x1(n)','psd of x2(n)');
grid on;
% end part (b) ------------------------------------------

%% Part (c)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];
w=filter(b,1,randn(64,1));

%Resolution
delw = 0.89*(2*pi/length(n));
disp('Resolution is:'),disp(delw);
disp('The values of x(n) to resolve frequencies in to two sinusoids:');
disp(length(n));

 x1=2.*cos(n.*w1+phi1);
 x2=2.*cos(n.*w2+phi2);
  x=x1+x2+w';
  
 pxx= periodogram(x);

 pxx_x1 = periodogram(x1);
 pxx_x2= periodogram(x2);

plot(10*log10(pxx));hold on
plot(10*log10(pxx_x1));hold on
plot(10*log10(pxx_x2));hold on

title('power spectrum');
xlabel('Frequency');
ylabel('psd in dB');
legend('psd of x(n)','psd of x1(n)','psd of x2(n)');
grid on;
 %%% The number of values of x(n) to solve the frequencies of the two
 %%% sinusoids are equal to length of phi.
%end part (c)--------------------------------------------------------

%% Part (d)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];
w=filter(b,1,randn(100,1));

 x1=2.*cos(n.*w1+phi1);
 x2=2.*cos(n.*w2+phi2);
  
 pxx= periodogram(w);
 
 pxx_x1 = periodogram(x1);
 pxx_x2= periodogram(x2);
 
plot(10*log10(pxx));hold on
plot(10*log10(pxx_x1));hold on
plot(10*log10(pxx_x2));hold on

title('power spectrum');
xlabel('Frequency');
ylabel('psd in dB');
legend('psd of w(n)','psd of x1(n)','psd of x2(n)');
grid on;
 %%% The number of values of x(n) to solve the frequencies of the two
 %%% sinusoids are equal to length of phi.
%end part (d)---------------------------------------------------------

%% Part (e)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];

 x1=2.*cos(n.*w1+phi1);
 x2=2.*cos(n.*w2+phi2);
 x_n1 =filter(x1,1,randn(256,1));
 x_n2 =filter(x2,1,randn(256,1));

 pxx_x1 = periodogram(x_n1);
 pxx_x2= periodogram(x_n2);
 figure(1)
subplot(2,1,1),plot(10*log10(pxx_x1));hold on
subplot(2,1,2),plot(10*log10(pxx_x2));hold on

title('power spectrum');
xlabel('Frequency');
ylabel('psd in dB');
%legend('psd of x1(n)','psd of x2(n)');
grid on;
 %%% The two spectrums are visible
%% For 20 different realizations
xn1 =filter(x1,1,randn(256,20));
xn2 =filter(x2,1,randn(256,20));
figure(2)

pxx_xn1 = periodogram(xn1);
pxx_xn2= periodogram(xn2);
 
subplot(2,1,1);plot(10*log10(pxx_xn1));
title('power spectrum of x1(n) for 20 realizations');
xlabel('Frequency');
ylabel('psd in dB');
grid on;
subplot(2,1,2);plot(10*log10(pxx_xn2));

title('power spectrum of x2(n) for 20 realizations');
xlabel('Frequency');
ylabel('psd in dB');
grid on;
%end part (e)----------------------------------------

%% Part (f)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];
K = 8;
N = 256;
L = N/K;

 x1=2.*cos(n.*w1+phi1);
 x2=2.*cos(n.*w2+phi2);
 x_n1 =filter(x1,1,randn(L,1));
 x_n2 =filter(x2,1,randn(L,1));

 pxx_x1 = bartlett(L);
 pxx_x2= bartlett(L);
 figure(1)
plot(10*log10(pxx_x1));hold on
plot(10*log10(pxx_x2));hold on

title('power spectrum');
xlabel('Frequency');
ylabel('psd in dB');
legend('psd of x1(n)','psd of x2(n)');
grid on;
 %%% The two spectrums are visible
%% For 20 different realizations
xn1 =filter(x1,1,randn(L,20));
xn2 =filter(x2,1,randn(L,20));
%end part(f)-------------------------------