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

00 reviews
0

Report Abuse

Description

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)-------------------------------

Power Spectrum

Power Spectrum

Power Spectrum

Author Profile

Mamush

Member since 4 years ago
    View Profile

    Contact Author

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

    Login to Write Your Review

    There are no reviews yet.