DSPRelated.com
   Code Snippets

PSD estimation using N-point DFT -scilab code

Senthilkumar ● December 28, 2011 ● Coded in Scilab 
//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')

Comparison of different power spectrum estimates

Senthilkumar ● December 28, 2011 ● Coded in Scilab 
//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

PSD Estimation [Blackmann-Tukey Method]-scilab code

Senthilkumar ● December 28, 2011 ● Coded in Scilab 
//Caption: PSD Estimation using Blackmann-Tukey Method
// of Given Discrete Sequence 
//x(n) =(0.9)^n, 0<= n<=20
clear;
clc;
close;
N =16; //8-point DFT
n = 0:20;
x1 = (0.9.^n) //given discrete sequence
x = x1([1:16]);
X = dft(x,-1); //16-point DFT of given discrete sequence
Pxx = (1/N)*(abs(X).^2); //Peridogram Estimate 
disp(x,'16 point input discrete sequence x(n)=')
disp(X,'DFT of x(n)is X(k)=')
disp(Pxx,'Peridogram of x(n) is Pxx(k/N)=')
figure(1)
a = gca();
a.data_bounds =[0,0;20,5];
plot2d3('gnn',[1:N],Pxx)
a.foreground = 5; 
a.font_color = 5;
a.font_style = 5;
title('Peridogram Estimate')
xlabel('Discrete Frequency Variable K ----->')
ylabel('Periodogram Pxx (k /N) ---->')
//Result
//16 point input discrete sequence x(n)=   
//column 1 to 9
//1.    0.9    0.81    0.729    0.6561    0.59049    0.531441    0.4782969    0.4304672  
//column 10 to 16
//0.3874205    0.3486784    0.3138106    0.2824295    0.2541866    0.2287679    0.2058911  
//DFT of x(n)is X(k)=   
//    8.1469798               
//    0.9337942 - 1.9085859i  
//    0.5514204 - 0.9651212i  
//    0.4763807 - 0.6042033i  
//    0.4501094 - 0.4050984i  
//    0.4383220 - 0.2710927i  
//    0.4324549 - 0.1681821i  
//    0.4296342 - 0.0807934i  
//    0.4287884 - 1.159D-15i  
//    0.4296342 + 0.0807934i  
//    0.4324549 + 0.1681821i  
//    0.4383220 + 0.2710927i  
//    0.4501094 + 0.4050984i  
//    0.4763807 + 0.6042033i  
//    0.5514204 + 0.9651212i  
//    0.9337942 + 1.9085859i  
// 
//Peridogram of x(n) is Pxx(k/N)=   
//    4.14833    
//    0.2821670  
//    0.0772202  
//    0.0370000  
//    0.0229190  
//    0.0166011  
//    0.0134564  
//    0.0119446  
//    0.0114912  
//    0.0119446  
//    0.0134564  
//    0.0166011  
//    0.0229190  
//    0.0370000  
//    0.0772202  
//    0.2821670

Subband coding speech signal-scilab code

Senthilkumar ● December 28, 2011●3 comments ● Coded in Scilab 
//Caption: write a program for subband coding of speech signal
clear;
clc;
[x,Fs,bits]=wavread("E:\4.wav");
n = length(x)
//Low Pass Filter of Length = 19, Wc = 0.5, Hamming Window
[wft_LPF,wfm_LPF,fr_LPF]=wfir('lp',18,[0.25,0],'hm',[0,0])
//High Pass Filter of Length = 19, Wc = 0.5, Hamming Window
[wft_HPF,wfm_HPF,fr_HPF]=wfir('hp',18,[0.25,0],'hm',[0,0])
//LPF output
Y_lpf = convol(x,wft_LPF)
//HPF output
Y_hpf = convol(x,wft_HPF)
//Downsampling by a factor of 2
Downsampling_XLPF = Y_lpf(1:2:length(Y_lpf));
Downsampling_XHPF = Y_hpf(1:2:length(Y_hpf));
figure(1)
subplot(3,1,1)
plot([1:n],x)
xtitle("ORIGINAL Speech SIGNAL");
subplot(3,1,2)
plot([1:length(Y_lpf)],Y_lpf)
xtitle("Low Pass Filtered Speech SIGNAL");
subplot(3,1,3)
plot([1:length(Y_hpf)],Y_hpf)
xtitle("High Pass Filtered Speech SIGNAL");
figure(2)
subplot(2,1,1)
plot([1:length(Y_lpf)],Y_lpf)
xtitle("Low Pass Filtered Speech SIGNAL");
subplot(2,1,2)
plot([1:length(Downsampling_XLPF)],Downsampling_XLPF,'r')
xtitle("LPF Output Downsampled by 2");
figure(3)
subplot(2,1,1)
plot([1:length(Y_hpf)],Y_hpf)
xtitle("High Pass Filtered Speech SIGNAL");
subplot(2,1,2)
plot([1:length(Downsampling_XHPF)],Downsampling_XHPF,'r')
xtitle("HPF Output Downsampled by 2");

Speech Noise cancellation - scilab code

Senthilkumar ● December 28, 2011●6 comments ● Coded in Scilab 
//Caption: Speech Noise Cancellation using LMS Adaptive Filter

clc;
//Reading a speech signal
[x,Fs,bits]=wavread("E:\4.wav");
order = 40;  // Adaptive filter order
x = x';
N = length(x);
t = 1:N;
//Plot the speech signal
figure(1)
subplot(2,1,1)
plot(t,x)
title('Noise free Speech Signal')
//Generation of noise signal
noise = 0.1*rand(1,length(x));
//Adding noise with speech signal
for i = 1:length(noise)
    primary(i)= x(i)+noise(i);
end
//Plot the noisy speech signal
subplot(2,1,2)
plot(t,primary)
title('primary = speech+noise (input 1)')
//Reference noise generation
for i = 1:length(noise)
    ref(i)= noise(i)+0.025*rand(10);
end
//Plot the reference noise
figure(2)
subplot(2,1,1)
plot(t,ref)
title('reference noise (input 2)')
//Adaptive filter coefficients initialized to zeros
w = zeros(order,1);
Average_Power = pow_1(x,N)
mu = 1/(10*order*Average_Power); //Adaptive filter step size
//Speech noise cancellation
for k = 1:110
    for i =1:N-order-1
        buffer = ref(i:i+order-1); //current order points of reference
        desired(i) = primary(i)-buffer'*w; // dot product the reference & filter
        w = w+(buffer.*mu*desired(i)); //update filter coefficients
    end
end
//Plot the Adaptive Filter output
subplot(2,1,2)
plot([1:length(desired)],desired)
title('Denoised Speech Signal at Adaptive Filter Output')
//Calculation of Mean Squarred Error between the original speech signal and
//Adaptive filter output
for i =1:N-order-1
    err(i) = x(i)-desired(i);
    square_error(i)= err(i)*err(i);
end
MSE = (sum(square_error))/(N-order-1);
MSE_dB = 20*log10(MSE);
//Playing the original speech signal
sound(x,Fs,16)
//Delay between playing sound signals
for i = 1:1000
    j = 1;
end
/////////////////////////////////
//Playing Noisy Speech Signal
sound(primary,Fs,16)
//Delay between playing sound signals
for i = 1:1000
    j = 1;
end
/////////////////////////////////
//Playing denoised speech signal (Adaptive Filter Output)
sound(desired,Fs,16)

speech signal-sampling rate conversion

Senthilkumar ● December 28, 2011●1 comment ● Coded in Scilab 
//Caption: write a program for decimation and interpolation of given
//Speech Signal [ Multirate Signal Processing]

clear;
clc;
[x,Fs,bits]=wavread("E:\4.wav");
n = length(x)
DECIMATION_OF_X = x(1:2:length(x));
INTERPOLATION_OF_X = zeros(1,2*length(x));
INTERPOLATION_OF_X(1:2:2*length(x)) = x;
subplot(3,1,1)
plot([1:n],x)
xtitle("ORIGINAL Speech SIGNAL");
subplot(3,1,2)
plot([1:ceil((n/2))],DECIMATION_OF_X)
xtitle("DECIMATION BY A FACTOR OF 2")
subplot(3,1,3)
plot([1:2*n],INTERPOLATION_OF_X)
xtitle("INTERPOLATION BY A FACTOR OF 2")

scilab program to speech signal information

Senthilkumar ● December 28, 2011●1 comment ● Coded in Scilab 
//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

sound play command in scilab

Senthilkumar ● December 28, 2011 ● Coded in Scilab 
//Caption: Reading a Speech Signal & 
//Play the signal
clear;
clc;
[y,Fs,bits_y]=wavread("E:\4.wav");
sound(y,Fs,16)

Plot the Low pass and High pass filtered speech signal

Senthilkumar ● December 28, 2011 ● Coded in Scilab 
//Caption: Program to Plot the frequency reponse of 
//[1]. Low Pass Filter 
//[2]. High Pass Filter

clear;
clc;
[x,Fs,bits]=wavread("E:\4.wav");
n = length(x)
//Low Pass Filter of Length = 19, Wc = 0.5, Hamming Window
[wft_LPF,wfm_LPF,fr_LPF]=wfir('lp',18,[0.25,0],'hm',[0,0])
//High Pass Filter of Length = 19, Wc = 0.5, Hamming Window
[wft_HPF,wfm_HPF,fr_HPF]=wfir('hp',18,[0.25,0],'hm',[0,0])

subplot(3,1,1)
plot([1:n],x)
xtitle("ORIGINAL Speech SIGNAL");
subplot(3,1,2)
a = gca();
plot(2*fr_LPF,wfm_LPF,'g')
poly1= a.children.children(1);
poly1.thickness =3;
xtitle("Frequency Response of Low Pass Filter N = 19, Wc =0.5, Hamming Window");
xgrid(1)
subplot(3,1,3)
b = gca();
plot(2*fr_HPF,wfm_HPF,'r')
poly1= b.children.children(1);
poly1.thickness =3;
xtitle("Frequency Response of High Pass Filter N= 19, Wc =0.5, Hamming Window");
xgrid(1)

Function used to calculate average power

Senthilkumar ● December 28, 2011 ● Coded in Scilab 
function[y]= pow_1(x,N)
    xold = 0.0;
    for n =1:N
        sumx = xold+(x(n)^2);
        xold = sumx;
    end
    y = sumx/N;
endfunction