DSPRelated.com

Direct Sequence Spread Spectrum (DS-BPSK)

Senthilkumar March 28, 20121 comment Coded in Scilab
function[st,mt]= DS_Spread_Spectrum(bt,ct_polar) 
//Caption:Direct Sequence Spread Coherent BPSK
//Generation of waveforms in DS/BPSK spread spectrum transmitter
//bt: Input Data Sequence (bipolar format)
//ct_polar: Spreading code (bipolar format)
Ft = 0:0.01:1; 
//bt = [1*ones(1,N) -1*ones(1,N)];
t = 0:length(bt)-1;
//ct_polar = [-1,-1,1,1,1,-1,1,-1,-1,1,1,1,-1,1];
mt = bt.*ct_polar;
Carrier = 2*sin(Ft*2*%pi);
st = [];
for i = 1:length(mt)
  st = [st mt(i)*Carrier];
end
//
figure
subplot(3,1,1)
a =gca();
a.x_location = "origin";
a.y_location = "origin";
a.data_bounds = [0,-2;20,2];
plot2d2(t,bt,5)
xlabel('                                                               t')
title('Data b(t)')
subplot(3,1,2)
a =gca();
a.x_location = "origin";
a.y_location = "origin";
a.data_bounds = [0,-2;20,2];
plot2d2(t,ct_polar,5)
xlabel('                                                                t')
title('Spreading code c(t)')
subplot(3,1,3)
a =gca();
a.x_location = "origin";
a.y_location = "origin";
a.data_bounds = [0,-2;20,2];
plot2d2(t,mt,5)
xlabel('                                                               t')
title('Product Signal m(t)')
//
figure
subplot(3,1,1)
a =gca();
a.x_location = "origin";
a.y_location = "origin";
a.data_bounds = [0,-2;20,2];
plot2d2(t,mt,5)
xlabel('                                                                t')
title('Product Signal m(t)')
subplot(3,1,2)
a =gca();
a.x_location = "origin";
a.y_location = "origin";
a.data_bounds = [0,-2;20,2];
plot(Carrier)
xlabel('                                                               t')
title('Carrier Signal')
subplot(3,1,3)
a =gca();
a.x_location = "origin";
a.y_location = "origin";
a.data_bounds = [0,-2;20,2];
plot(st)
xlabel('                                                               t')
title('DS/BPSK signal s(t)')
endfunction
//Result
//->bt = [1,1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,-1]
// bt  =
// 
//    1.    1.    1.    1.    1.    1.    1.  - 1.  - 1.  - 1.  - 1.  - 1.  - 1.  - 1.  
// 
//-->ct_polar = [-1,-1,1,1,1,-1,1,-1,-1,1,1,1,-1,1]
// ct_polar  =
// 
//  - 1.  - 1.    1.    1.    1.  - 1.    1.  - 1.  - 1.    1.    1.    1.  - 1.    1.  
// 
//-->exec('C:\Users\SENTHILKUMAR\Desktop\Communication_Toolbox\Digital_Communication\DS_Spread_Spectrum.sci', -1)
// 
//-->[st,mt]= DS_Spread_Spectrum(bt,ct_polar)
// mt  =
// 
//  - 1.  - 1.    1.    1.    1.  - 1.    1.    1.    1.  - 1.  - 1.  - 1.    1.  - 1.

Convolutional coding using transform domain approach

Senthilkumar March 28, 2012 Coded in Scilab
function [x1D,x2D]= ConvolutionCode_TransDomain()         
//Convolutional code - Transform domain approach
//g1D =  generator polynomial 1
//g2D =  generator polynomial 2
//x1D = top output sequence polynomial
//x2D = bottom output sequence polynomial 
D = poly(0,'D');
g1D = input('Enter the generator polynomial 1=') //generator polynomial 1
g2D = input('Enter the generator polynomial 2=') //generator polynomial 2
mD = input('Enter the message sequence')//message sequence polynomial representation
x1D = g1D*mD; //top output polynomial
x2D = g2D*mD; //bottom output polynomial
x1 = coeff(x1D);
x2 = coeff(x2D);
disp(modulo(x1,2),'top output sequence')
disp(modulo(x2,2),'bottom output sequence')
endfunction
//Result
//Enter the generator polynomial 1 =  1+D+D^2
//Enter the generator polynomial 2 =  1+D^2;
//Enter the message sequence = 1+0+0+D^3+D^4;
//top output sequence   
// 
//    1.    1.    1.    1.    0.    0.    1.   
//bottom output sequence   
// 
//    1.    0.    1.    1.    1.    1.    1.  
//x2D  =
// 
//         2   3   4   5   6  
//    1 + D + D + D + D + D   
//x1D  =
// 
//             2   3    4    5   6  
//    1 + D + D + D + 2D + 2D + D

Constellation Diagram for Binary QPSK

Senthilkumar March 28, 20122 comments Coded in Scilab
function[y] = Constellation_QPSK()
M =4;
i = 1:M;
y = cos((2*i-1)*%pi/4)-sin((2*i-1)*%pi/4)*%i;
annot = dec2bin([0:M-1],log2(M));
disp(y,'coordinates of message points')
disp(annot,'dibits value')
figure;
a =gca();
a.data_bounds = [-1,-1;1,1];
a.x_location = "origin";
a.y_location = "origin";
plot2d(real(y(1)),imag(y(1)),-2)
plot2d(real(y(2)),imag(y(2)),-4)
plot2d(real(y(3)),imag(y(3)),-5)
plot2d(real(y(4)),imag(y(4)),-9)
xlabel('                                             In-Phase');
ylabel('                                             Quadrature');
title('Constellation for QPSK')
legend(['message point 1 (dibit 10)';'message point 2 (dibit 00)';'message point 3 (dibit 01)';'message point 4 (dibit 11)'],5)
endfunction
//Result
//coordinates of message points   
//  0.7071068 - 0.7071068i  - 0.7071068 - 0.7071068i  - 0.7071068 + 0.7071068i    0.7071068 + 0.7071068i  
// 
// dibits value   
// 
//!00  01  10  11  !

Constellation diagram for Binary PSK

Senthilkumar March 28, 2012 Coded in Scilab
function[y]= Constellation_BPSK()
M =2;
i = 1:M;
y = cos(2*%pi+(i-1)*%pi);
annot = dec2bin([length(y)-1:-1:0],log2(M));
disp(y,'coordinates of message points')
disp(annot,'Message points')
figure;
a =gca();
a.data_bounds = [-2,-2;2,2];
a.x_location = "origin";
a.y_location = "origin";
plot2d(real(y(1)),imag(y(1)),-9)
plot2d(real(y(2)),imag(y(2)),-5)
xlabel('                                                                      In-Phase');
ylabel('                                                                      Quadrature');
title('Constellation for BPSK')
legend(['message point 1 (binary 1)';'message point 2 (binary 0)'],5)
endfunction
//Result
//coordinates of message points   
// 
//    1.  - 1.  
// 
// Message points   
// 
//!1  0  !

Constellation for Binary FSK

Senthilkumar March 28, 2012 Coded in Scilab
function[y]= Constellation_BFSK()
M =2;
y = [1,0;0,1];
annot = dec2bin([M-1:-1:0],log2(M));
disp(y,'coordinates of message points')
disp(annot,'Message points')
figure;
a =gca();
a.data_bounds = [-2,-2;2,2];
a.x_location = "origin";
a.y_location = "origin";
plot2d(y(1,1),y(1,2),-9)
plot2d(y(2,1),y(2,2),-5)
xlabel('                                                                      In-Phase');
ylabel('                                                                    Quadrature');
title('Constellation for BFSK')
legend(['message point 1 (binary 1)';'message point 2 (binary 0)'],5)
endfunction
//Result
//coordinates of message points   
// 
//    1.    0.  
//    0.    1.  
// 
// Message points   
// 
//!1  0  !

PSD estimation window based-scilab code

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

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, 20113 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");