//Caption: Evaluating power spectrum of a discrete sequence
//Using N-point DFT
clear all;
clc;
close;
N =16; //Number of samples in given sequence
n =0:N-1;
delta_f = [0.06,0.01];//frequency separation
x1 = sin(2*%pi*0.315*n)+cos(2*%pi*(0.315+delta_f(1))*n);
x2 = sin(2*%pi*0.315*n)+cos(2*%pi*(0.315+delta_f(2))*n);
L = [8,16,32,128];
k1 = 0:L(1)-1;
k2 = 0:L(2)-1;
k3 = 0:L(3)-1;
k4 = 0:L(4)-1;
fk1 = k1./L(1);
fk2 = k2./L(2);
fk3 = k3./L(3);
fk4 = k4./L(4);
for i =1:length(fk1)
Pxx1_fk1(i) = 0;
Pxx2_fk1(i) = 0;
for m = 1:N
Pxx1_fk1(i)=Pxx1_fk1(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk1(i));
Pxx2_fk1(i)=Pxx1_fk1(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk1(i));
end
Pxx1_fk1(i) = (Pxx1_fk1(i)^2)/N;
Pxx2_fk1(i) = (Pxx2_fk1(i)^2)/N;
end
for i =1:length(fk2)
Pxx1_fk2(i) = 0;
Pxx2_fk2(i) = 0;
for m = 1:N
Pxx1_fk2(i)=Pxx1_fk2(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk2(i));
Pxx2_fk2(i)=Pxx1_fk2(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk2(i));
end
Pxx1_fk2(i) = (Pxx1_fk2(i)^2)/N;
Pxx2_fk2(i) = (Pxx1_fk2(i)^2)/N;
end
for i =1:length(fk3)
Pxx1_fk3(i) = 0;
Pxx2_fk3(i) = 0;
for m = 1:N
Pxx1_fk3(i) =Pxx1_fk3(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk3(i));
Pxx2_fk3(i) =Pxx1_fk3(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk3(i));
end
Pxx1_fk3(i) = (Pxx1_fk3(i)^2)/N;
Pxx2_fk3(i) = (Pxx1_fk3(i)^2)/N;
end
for i =1:length(fk4)
Pxx1_fk4(i) = 0;
Pxx2_fk4(i) = 0;
for m = 1:N
Pxx1_fk4(i) =Pxx1_fk4(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk4(i));
Pxx2_fk4(i) =Pxx1_fk4(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk4(i));
end
Pxx1_fk4(i) = (Pxx1_fk4(i)^2)/N;
Pxx2_fk4(i) = (Pxx1_fk4(i)^2)/N;
end
figure(1)
subplot(2,2,1)
plot2d3('gnn',k1,abs(Pxx1_fk1))
xtitle('8 point DFT')
subplot(2,2,2)
plot2d3('gnn',k2,abs(Pxx1_fk2))
xtitle('16 point DFT')
subplot(2,2,3)
plot2d3('gnn',k3,abs(Pxx1_fk3))
xtitle('32 point DFT')
subplot(2,2,4)
plot2d3('gnn',k4,abs(Pxx1_fk4))
xtitle('128 point DFT')
figure(2)
subplot(2,2,1)
plot2d3('gnn',k1,abs(Pxx2_fk1))
xtitle('8 point DFT')
subplot(2,2,2)
plot2d3('gnn',k2,abs(Pxx2_fk2))
xtitle('16 point DFT')
subplot(2,2,3)
plot2d3('gnn',k3,abs(Pxx2_fk3))
xtitle('32 point DFT')
subplot(2,2,4)
plot2d3('gnn',k4,abs(Pxx2_fk4))
xtitle('128 point DFT')
//Caption:Determination of spectrum of a signal
//With maximum normalized frequency f = 0.1
//using Rectangular window and Blackmann window
clear all;
close;
clc;
N = 61;
cfreq = [0.1 0];
[wft,wfm,fr]=wfir('lp',N,cfreq,'re',0);
disp(wft,'Time domain filter coefficients hd(n)=');
disp(wfm,'Frequency domain filter values Hd(w)=');
WFM_dB = 20*log10(wfm);//Frequency response in dB
for n = 1:N
h_balckmann(n)=0.42-0.5*cos(2*%pi*n/(N-1))+0.08*cos(4*%pi*n/(N-1));
end
wft_blmn = wft'.*h_balckmann;
disp(wft_blmn,'Blackmann window based Filter output h(n)=')
wfm_blmn = frmag(wft_blmn,length(fr));
WFM_blmn_dB =20*log10(wfm_blmn);
subplot(2,1,1)
plot2d(fr,WFM_dB)
xgrid(1)
xtitle('Power Spectrum with Rectangular window Filtered M = 61','Frequency in cycles per samples f','Energy density in dB')
subplot(2,1,2)
plot2d(fr,WFM_blmn_dB)
xgrid(1)
xtitle('Power Spectrum with Blackmann window Filtered M = 61','Frequency in cycles per samples f','Energy density in dB')
//Autocorrelation of a given Input Sequence
//Finding out the period of the signal using autocorrelation technique
clear;
clc;
close;
x = input('Enter the given discrete time sequence');
L = length(x);
h = zeros(1,L);
for i = 1:L
h(L-i+1) = x(i);
end
N = 2*L-1;
Rxx = zeros(1,N);
for i = L+1:N
h(i) = 0;
end
for i = L+1:N
x(i) = 0;
end
for n = 1:N
for k = 1:N
if(n >= k)
Rxx(n) = Rxx(n)+x(n-k+1)*h(k);
end
end
end
disp('Auto Correlation Result is')
Rxx
disp('Center Value is the Maximum of autocorrelation result')
[m,n] = max(Rxx)
disp('Period of the given signal using Auto Correlation Sequence')
n
//Caption: Reading a Speech Signal &
//[1]. Displaying its sampling rate
//[2]. Number of bits used per speech sample
//[3]. Total Time duration of the speech signal in seconds
clear;
clc;
[y,Fs,bits]=wavread("E:\4.wav");
a = gca();
plot2d(y);
a.x_location = 'origin';
title('Speech signal with Sampling Rate = 8 KHz, No. of Samples = 8360')
disp(Fs,'Sampling Rate in Hz Fs = ');
disp(bits,'Number of bits used per speech sample b =');
N = length(y);
T = N/Fs;
disp(N,'Total Number of Samples N =')
disp(T,'Duration of speech signal in seconds T=')
//Result
//Sampling Rate in Hz Fs =
// 8000.
//Number of bits used per speech sample b =
// 16.
//Total Number of Samples N =
// 8360.
//Duration of speech signal in seconds T=
// 1.045
//Continuous Time Fourier Series Spectral Coefficients of
//a periodic sine signal x(t) = sin(Wot)
clear;
close;
clc;
t = 0:0.01:1;
T = 1;
Wo = 2*%pi/T;
xt = sin(Wo*t);
for k =0:5
C(k+1,:) = exp(-sqrt(-1)*Wo*t.*k);
a(k+1) = xt*C(k+1,:)'/length(t); //fourier series is done
if(abs(a(k+1))<=0.01)
a(k+1)=0;
end
end
a =a';
ak = [-a($:-1:1),a(2:$)];
disp(ak,'Continuous Time Fourier Series Coefficients are:')
//Result
//Continuous Time Fourier Series Coefficients are:
// column 1 to 11
// 0 0 0 0 -1.010D-18+0.4950495i 0 1.010D-18-0.4950495i 0 0 0 0
//Caption:Determination of spectrum of a signal
//With maximum normalized frequency f = 0.1
//using Rectangular window and Blackmann window
clear all;
close;
clc;
N = 61;
cfreq = [0.1 0];
[wft,wfm,fr]=wfir('lp',N,cfreq,'re',0);
disp(wft,'Time domain filter coefficients hd(n)=');
disp(wfm,'Frequency domain filter values Hd(w)=');
WFM_dB = 20*log10(wfm);//Frequency response in dB
for n = 1:N
h_balckmann(n)=0.42-0.5*cos(2*%pi*n/(N-1))+0.08*cos(4*%pi*n/(N-1));
end
wft_blmn = wft'.*h_balckmann;
disp(wft_blmn,'Blackmann window based Filter output h(n)=')
wfm_blmn = frmag(wft_blmn,length(fr));
WFM_blmn_dB =20*log10(wfm_blmn);
subplot(2,1,1)
plot2d(fr,WFM_dB)
xgrid(1)
xtitle('Power Spectrum with Rectangular window Filtered M = 61','Frequency in cycles per samples f','Energy density in dB')
subplot(2,1,2)
plot2d(fr,WFM_blmn_dB)
xgrid(1)
xtitle('Power Spectrum with Blackmann window Filtered M = 61','Frequency in cycles per samples f','Energy density in dB')
//Caption: Evaluating power spectrum of a discrete sequence
//Using N-point DFT
clear all;
clc;
close;
N =16; //Number of samples in given sequence
n =0:N-1;
delta_f = [0.06,0.01];//frequency separation
x1 = sin(2*%pi*0.315*n)+cos(2*%pi*(0.315+delta_f(1))*n);
x2 = sin(2*%pi*0.315*n)+cos(2*%pi*(0.315+delta_f(2))*n);
L = [8,16,32,128];
k1 = 0:L(1)-1;
k2 = 0:L(2)-1;
k3 = 0:L(3)-1;
k4 = 0:L(4)-1;
fk1 = k1./L(1);
fk2 = k2./L(2);
fk3 = k3./L(3);
fk4 = k4./L(4);
for i =1:length(fk1)
Pxx1_fk1(i) = 0;
Pxx2_fk1(i) = 0;
for m = 1:N
Pxx1_fk1(i)=Pxx1_fk1(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk1(i));
Pxx2_fk1(i)=Pxx1_fk1(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk1(i));
end
Pxx1_fk1(i) = (Pxx1_fk1(i)^2)/N;
Pxx2_fk1(i) = (Pxx2_fk1(i)^2)/N;
end
for i =1:length(fk2)
Pxx1_fk2(i) = 0;
Pxx2_fk2(i) = 0;
for m = 1:N
Pxx1_fk2(i)=Pxx1_fk2(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk2(i));
Pxx2_fk2(i)=Pxx1_fk2(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk2(i));
end
Pxx1_fk2(i) = (Pxx1_fk2(i)^2)/N;
Pxx2_fk2(i) = (Pxx1_fk2(i)^2)/N;
end
for i =1:length(fk3)
Pxx1_fk3(i) = 0;
Pxx2_fk3(i) = 0;
for m = 1:N
Pxx1_fk3(i) =Pxx1_fk3(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk3(i));
Pxx2_fk3(i) =Pxx1_fk3(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk3(i));
end
Pxx1_fk3(i) = (Pxx1_fk3(i)^2)/N;
Pxx2_fk3(i) = (Pxx1_fk3(i)^2)/N;
end
for i =1:length(fk4)
Pxx1_fk4(i) = 0;
Pxx2_fk4(i) = 0;
for m = 1:N
Pxx1_fk4(i) =Pxx1_fk4(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk4(i));
Pxx2_fk4(i) =Pxx1_fk4(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk4(i));
end
Pxx1_fk4(i) = (Pxx1_fk4(i)^2)/N;
Pxx2_fk4(i) = (Pxx1_fk4(i)^2)/N;
end
figure(1)
subplot(2,2,1)
plot2d3('gnn',k1,abs(Pxx1_fk1))
xtitle('8 point DFT')
subplot(2,2,2)
plot2d3('gnn',k2,abs(Pxx1_fk2))
xtitle('16 point DFT')
subplot(2,2,3)
plot2d3('gnn',k3,abs(Pxx1_fk3))
xtitle('32 point DFT')
subplot(2,2,4)
plot2d3('gnn',k4,abs(Pxx1_fk4))
xtitle('128 point DFT')
figure(2)
subplot(2,2,1)
plot2d3('gnn',k1,abs(Pxx2_fk1))
xtitle('8 point DFT')
subplot(2,2,2)
plot2d3('gnn',k2,abs(Pxx2_fk2))
xtitle('16 point DFT')
subplot(2,2,3)
plot2d3('gnn',k3,abs(Pxx2_fk3))
xtitle('32 point DFT')
subplot(2,2,4)
plot2d3('gnn',k4,abs(Pxx2_fk4))
xtitle('128 point DFT')
//Caption: Determination of Frequency Resolution of the
//[1]. Barlett [2]. Welch [3].Blackman Tukey Power Spectrum Estimation
//Methods
clear;
clc;
close;
Q = 10; //Quality factor
N = 1000; //Length of the sample sequence
//Bartlett Method
F_Bartlett = Q/(1.11*N);
disp(F_Bartlett,'Frequency Resolution of Bartlett Power Spectrum Estimation')
//Welch Method
F_Welch = Q/(1.39*N);
disp(F_Welch,'Frequency Resolution of Welch Power Spectrum Estimation')
//Blackmann-Tukey Method
F_Blackmann_Tukey = Q/(2.34*N);
disp(F_Blackmann_Tukey,'Frequency Resolution of Blackmann Tukey Power Spectrum Estimation')
//Result
//Frequency Resolution of Bartlett Power Spectrum Estimation
// 0.0090090
//Frequency Resolution of Welch Power Spectrum Estimation
// 0.0071942
//Frequency Resolution of Blackmann Tukey Power Spectrum Estimation
// 0.0042735
