Code Snippets Submitted by rsenthil_signal
Constellation Diagram for Binary QPSK
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 !Discrete Plot in scilab
//This Program Illustrates the discrete plot in scilab
//using plot2d3 function
clear;
clc;
close;
a =1.5;
n =1:10;
x = (a)^n;
a=gca();
a.thickness = 2;
plot2d3('gnn',n,x)
xtitle('Graphical Representation of Exponentially Increasing Signal','n','x[n]');Speech Noise cancellation - scilab code
//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)Analog IIR Butterworth Filter - Scilab
// To Design an Analog Butterworth Filter
//For the given cutoff frequency and filter order
//Wc = 500 Hz
omegap = 500; //pass band edge frequency
omegas = 1000;//stop band edge frequency
delta1_in_dB = -3;//PassBand Ripple in dB
delta2_in_dB = -40;//StopBand Ripple in dB
delta1 = 10^(delta1_in_dB/20)
delta2 = 10^(delta2_in_dB/20)
//Caculation of filter order
N = log10((1/(delta2^2))-1)/(2*log10(omegas/omegap))
N = ceil(N) //Rounding off nearest integer
omegac = omegap;
h=buttmag(N,omegac,1:1000);//Analog Butterworth filter magnitude response
mag=20*log10(h);//Magntitude Response in dB
plot2d((1:1000),mag,[0,-180,1000,20]);
a=gca();
a.thickness = 3;
a.foreground = 1;
a.font_style = 9;
xgrid(5)
xtitle('Magnitude Response of Butterworth LPF Filter Cutoff frequency = 500 Hz','Analog frequency in Hz--->','Magnitude in dB -->');Decimation or Downsampling in scilab
//Multirate Signal Processing in scilab
//Downsampling a sinusoidal signal by a factor of 2
clear;
clc;
n = 0:%pi/200:2*%pi;
x = sin(%pi*n); //original signal
downsampling_x = x(1:2:length(x)); //downsampled by a factor of 2
subplot(2,1,1)
plot(1:length(x),x);
xtitle('original singal')
subplot(2,1,2)
plot(1:length(downsampling_x),downsampling_x);
xtitle('Downsampled Signal by a factor of 2');Duobinary Encoder and Decoder
function [c,b_r]=Duobinary_EncDec(b)
// Precoded Duobinary coder and decoder
//b = input binary sequence:precoder input
//c = duobinary coder output
//b_r = duobinary decoder output
a(1) = xor(1,b(1));
if(a(1)==1)
a_volts(1) = 1;
end
for k =2:length(b)
a(k) = xor(a(k-1),b(k));
if(a(k)==1)
a_volts(k)=1;
else
a_volts(k)=-1;
end
end
a = a';
a_volts = a_volts';
disp(a,'Precoder output in binary form:')
disp(a_volts,'Precoder output in volts:')
//Duobinary coder output in volts
c(1) = 1+ a_volts(1);
for k =2:length(a)
c(k) = a_volts(k-1)+a_volts(k);
end
c = c';
disp(c,'Duobinary coder output in volts:')
//Duobinary decoder output by applying decision rule
for k =1:length(c)
if(abs(c(k))>1)
b_r(k) = 0;
else
b_r(k) = 1;
end
end
b_r = b_r';
disp(b_r,'Recovered original sequence at detector oupupt:')
endfunction
//Result
//Precoder output in binary form:
//
// 1. 1. 0. 0. 1. 0. 0.
//
// Precoder output in volts:
//
// 1. 1. - 1. - 1. 1. - 1. - 1.
//
// Duobinary coder output in volts:
//
// 2. 2. 0. - 2. 0. 0. - 2.
//
// Recovered original sequence at detector oupupt:
//
// 0. 0. 1. 0. 1. 1. 0.NLMS Adaptive filter
%Normalized Least Mean Square Adaptive Filter
%for Echo Cancellation
function [e,h,t] = adapt_filt_nlms(x,B,h,delta,l,l1)
%x = the signal from the speaker
%B = signal through echo path
%h = impulse response of adaptive filter
%l = length of the signal x
%l1 = length of the adaptive filter
%delta = step size
%e = error
%h = adaptive filter impulse response
tic;
for k = 1:150
for n= 1:l
xcap = x((n+l1-1):-1:(n+l1-1)-l1+1);
yout(n) = h*xcap';
e(n) = B(n) - yout(n);
xnorm = 0.000001+(xcap*xcap');
h = h+((delta*e(n))*(xcap/xnorm));
end
eold = 0.0;
for i = 1:l
mesquare = eold+(e(i)^2);
eold = mesquare;
end
if mesquare <=0.0001
break;
end
end
t = toc; %to get the time elapsedu -Law and Inverse u-Law
function [Cx,Xmax] = mulaw(x,mu)
//Non-linear Quantization
//mulaw: mulaw nonlinear quantization
//x = input vector
//Cx = mulaw compressor output
//Xmax = maximum of input vector x
Xmax = max(abs(x));
if(log(1+mu)~=0)
Cx = (log(1+mu*abs(x/Xmax))./log(1+mu));
else
Cx = x/Xmax;
end
Cx = Cx/Xmax; //normalization of output vector
endfunction
////////////////////////////////////////////////////
function x = invmulaw(y,mu)
//Non-linear Quantization
//invmulaw: inverse mulaw nonlinear quantization
//x = output vector
//y = input vector (using mulaw nonlinear comression)
x = (((1+mu).^(abs(y))-1)./mu)
endfunctionCircular convolution
//Performing Circular COnvolution
//Using DFT
clear all;
clc;
close;
L = 4; //Length of the Sequence
N = 4; // N -point DFT
x1 = [2,1,2,1];
x2 = [1,2,3,4];
//Computing DFT
X1 = dft(x1,-1)
X2 = dft(x2,-1)
//Multiplication of 2 DFTs
X3 = X1.*X2
//Circular Convolution Result
x3 =abs(dft(X3,1))
disp(x3,'x3=')Image Denoising - Bayesshrink
%Using BAYESSHRINK
%subband dependent threshold
%Soft threshold
function [soft_X1,SOFT_PSNR] = Bayes_soft(X,Y)
%One -level decomposition
[CA,CH,CV,CD] = dwt2(Y,'haar');
%Call the function to calculate the threshold
CH =CH/1.2;
CV = CV/1.2;
CD = CD/1.2;
T1_CH = bayeshrink(CD,CH);
T1_CV = bayeshrink(CD,CV);
T1_CD = bayeshrink(CD,CD);
% Call the function to perfom soft shrinkage
de_CH = soft1(CH,T1_CH);
de_CV = soft1(CV,T1_CV);
de_CD = soft1(CD,T1_CD);
%Two -level decomposition
[CA1,CH1,CV1,CD1] = dwt2(CA,'haar');
CH1 = CH1/1.2;
CV1 = CV1/1.2;
CD1 = CD1/1.2;
%Call the function to calculate the threshold
T1_CH1 = bayeshrink(CD1,CH1);
T1_CV1 = bayeshrink(CD1,CV1);
T1_CD1 = bayeshrink(CD1,CD1);
% Call the function to perfom soft shrinkage
de_CH1 = soft1(CH1,T1_CH1);
de_CV1 = soft1(CV1,T1_CV1);
de_CD1 = soft1(CD1,T1_CD1);
%Reconstruction for soft shrinkage
X2 = idwt2(CA1,de_CH1,de_CV1,de_CD1,'haar');
X1 = idwt2(X2,de_CH,de_CV,de_CD,'haar');
SOFT_PSNR = PSNR(X,X1);
soft_X1 = uint8(X1);Power spectra of MPSK
function[SB_PSK] = PowerSpectra_MPSK()
rb = input('Enter the bit rate=');
Eb = input('Enter the energy of the bit=');
f = 0:1/100:rb;
Tb = 1/rb; //Bit duration
M = [2,4,8];
for j = 1:length(M)
for i= 1:length(f)
SB_PSK(j,i)=2*Eb*(sinc_new(f(i)*Tb*log2(M(j)))^2)*log2(M(j));
end
end
a=gca();
plot2d(f*Tb,SB_PSK(1,:)/(2*Eb))
plot2d(f*Tb,SB_PSK(2,:)/(2*Eb),2)
plot2d(f*Tb,SB_PSK(3,:)/(2*Eb),5)
xlabel('Normalized Frequency ---->')
ylabel('Normalized Power Spectral Density--->')
title('Power Spectra of M-ary signals for M =2,4,8')
legend(['M=2','M=4','M=8'])
xgrid(1)
endfunction
//Result
//Enter the bit rate in bits per second:2
//Enter the Energy of bit:1Power spectra of MSK & QPSk
function [SB_MSK,SB_QPSK]= PowerSpectra_MSK_QPSK()
//Comparison of QPSK and MSK Power Spectrums
rb = input('Enter the bit rate in bits per second:');
Eb = input('Enter the Energy of bit:');
f = 0:1/(100*rb):(4/rb);
Tb = 1/rb; //bit duration in seconds
for i = 1:length(f)
if(f(i)==0.5)
SB_MSK(i) = 4*Eb*f(i);
else
SB_MSK(i) = (32*Eb/(%pi^2))*(cos(2*%pi*Tb*f(i))/((4*Tb*f(i))^2-1))^2;
end
SB_QPSK(i)= 4*Eb*sinc_new((2*Tb*f(i)))^2;
end
a = gca();
plot(f*Tb,SB_MSK/(4*Eb));
plot(f*Tb,SB_QPSK/(4*Eb));
poly1= a.children(1).children(1);
poly1.foreground = 3;
xlabel('Normalized Frequency ---->')
ylabel('Normalized Power Spectral Density--->')
title('QPSK Vs MSK Power Spectra Comparison')
legend(['Minimum Shift Keying','QPSK'])
xgrid(1)
endfunction
//Result
//Enter the bit rate in bits per second:2
//Enter the Energy of bit:1Power Spectrum of Discrete PAM signals
function [Sxxf_NRZ_P,Sxxf_NRZ_BP,Sxxf_NRZ_UP,Sxxf_Manch]=PowerSpectra_PAM()
a = input('Enter the Amplitude value:');
fb = input('Enter the bit rate:');
Tb = 1/fb; //bit duration
f = 0:1/(100*Tb):2/Tb;
for i = 1:length(f)
Sxxf_NRZ_P(i) = (a^2)*Tb*(sinc_new(f(i)*Tb)^2);
Sxxf_NRZ_BP(i) = (a^2)*Tb*((sinc_new(f(i)*Tb))^2)*((sin(%pi*f(i)*Tb))^2);
if (i==1)
Sxxf_NRZ_UP(i) = (a^2)*(Tb/4)*((sinc_new(f(i)*Tb))^2)+(a^2)/4;
else
Sxxf_NRZ_UP(i) = (a^2)*(Tb/4)*((sinc_new(f(i)*Tb))^2);
end
Sxxf_Manch(i) = (a^2)*Tb*(sinc_new(f(i)*Tb/2)^2)*(sin(%pi*f(i)*Tb/2)^2);
end
//Plotting
a = gca();
plot2d(f,Sxxf_NRZ_P)
poly1= a.children(1).children(1);
poly1.thickness = 2; // the tickness of a curve.
plot2d(f,Sxxf_NRZ_BP,2)
poly1= a.children(1).children(1);
poly1.thickness = 2; // the tickness of a curve.
plot2d(f,Sxxf_NRZ_UP,5)
poly1= a.children(1).children(1);
poly1.thickness = 2; // the tickness of a curve.
plot2d(f,Sxxf_Manch,9)
poly1= a.children(1).children(1);
poly1.thickness = 2; // the tickness of a curve.
xlabel('f*Tb------->')
ylabel('Sxx(f)------->')
title('Power Spectral Densities of Different Line Codinig Techniques')
xgrid(1)
legend(['NRZ Polar Format','NRZ Bipolar format','NRZ Unipolar format','Manchester format']);
endfunction
//Result
//Enter the Amplitude value:1
//Enter the bit rate:1Raised Cosine Spectrum
function [p]= RaisedCosineSpectrum()
//Practical Solution for Intersymbol Interference
//Raised Cosine Spectrum
rb = input('Enter the bit rate:');
Tb =1/rb;
t =-3:1/100:3;
Bo = rb/2;
Alpha =0; //roll-off factor Intialized to zero
x =t/Tb;
for j =1:3
for i =1:length(t)
if((j==3)&((t(i)==0.5)|(t(i)==-0.5)))
p(j,i) = sinc_new(2*Bo*t(i));
else
num = sinc_new(2*Bo*t(i))*cos(2*%pi*Alpha*Bo*t(i));
den = 1-16*(Alpha^2)*(Bo^2)*(t(i)^2)+0.01;
p(j,i)= num/den;
end
end
Alpha = Alpha+0.5;
end
a =gca();
plot2d(t,p(1,:))
plot2d(t,p(2,:))
poly1= a.children(1).children(1);
poly1.foreground=2;
plot2d(t,p(3,:))
poly2= a.children(1).children(1);
po1y2.foreground=4;
poly2.line_style = 3;
xlabel('t/Tb------>');
ylabel('p(t)------->');
title('RAISED COSINE SPECTRUM - Practical Solution for ISI')
legend(['ROlloff Factor =0','ROlloff Factor =0.5','ROlloff Factor =1'])
xgrid(1)
endfunction
//Result
//Enter the bit rate:1Modified Duobinary Signaling-Spectrum
function[Amplitude_Response,Phase_Response]=Modified_Duobinary_Signaling()
//Modified Duobinary Signaling Scheme
//Magnitude and Phase Response
rb = input('Enter the bit rate=');
Tb =1/rb; //Bit duration
f = -rb/2:1/100:rb/2;
Amplitude_Response = abs(2*sin(2*%pi*f.*Tb));
Phase_Response = -(2*%pi*f.*Tb);
subplot(2,1,1)
a=gca();
a.x_location ="origin";
a.y_location ="origin";
plot2d(f,Amplitude_Response,2)
poly1= a.children(1).children(1);
poly1.thickness = 2; // the tickness of a curve.
xlabel('Frequency f---->')
ylabel('|H(f)| ----->')
title('Amplitude Repsonse of Modified Duobinary Singaling')
xgrid(1)
subplot(2,1,2)
a=gca();
a.x_location ="origin";
a.y_location ="origin";
plot2d(f,Phase_Response,5)
poly1= a.children(1).children(1);
poly1.thickness = 2; // the tickness of a curve.
xlabel(' Frequency f---->')
ylabel(' <H(f) ----->')
title('Phase Repsonse of Modified Duobinary Singaling')
xgrid(1)
//Result
//Enter the bit rate=8
endfunction
//Result
//Enter the bit rate=8Duobinary Encoder and Decoder
function [c,b_r]=Duobinary_EncDec(b)
// Precoded Duobinary coder and decoder
//b = input binary sequence:precoder input
//c = duobinary coder output
//b_r = duobinary decoder output
a(1) = xor(1,b(1));
if(a(1)==1)
a_volts(1) = 1;
end
for k =2:length(b)
a(k) = xor(a(k-1),b(k));
if(a(k)==1)
a_volts(k)=1;
else
a_volts(k)=-1;
end
end
a = a';
a_volts = a_volts';
disp(a,'Precoder output in binary form:')
disp(a_volts,'Precoder output in volts:')
//Duobinary coder output in volts
c(1) = 1+ a_volts(1);
for k =2:length(a)
c(k) = a_volts(k-1)+a_volts(k);
end
c = c';
disp(c,'Duobinary coder output in volts:')
//Duobinary decoder output by applying decision rule
for k =1:length(c)
if(abs(c(k))>1)
b_r(k) = 0;
else
b_r(k) = 1;
end
end
b_r = b_r';
disp(b_r,'Recovered original sequence at detector oupupt:')
endfunction
//Result
//Precoder output in binary form:
//
// 1. 1. 0. 0. 1. 0. 0.
//
// Precoder output in volts:
//
// 1. 1. - 1. - 1. 1. - 1. - 1.
//
// Duobinary coder output in volts:
//
// 2. 2. 0. - 2. 0. 0. - 2.
//
// Recovered original sequence at detector oupupt:
//
// 0. 0. 1. 0. 1. 1. 0.Duobianry Signaling-Amplitude & Phase Response
function [Amplitude_Response,Phase_Response]= Duobinary_Signaling()
//Duobinary Signaling Scheme
//Magnitude and Phase Response
rb = input('Enter the bit rate=');
Tb =1/rb; //Bit duration
f = -rb/2:1/100:rb/2;
Amplitude_Response = abs(2*cos(%pi*f.*Tb));
Phase_Response = -(%pi*f.*Tb);
subplot(2,1,1)
a=gca();
a.x_location ="origin";
a.y_location ="origin";
plot(f,Amplitude_Response)
xlabel('Frequency f---->')
ylabel('|H(f)| ----->')
title('Amplitude Repsonse of Duobinary Singaling')
subplot(2,1,2)
a=gca();
a.x_location ="origin";
a.y_location ="origin";
plot(f,Phase_Response)
xlabel(' Frequency f---->')
ylabel(' <H(f) ----->')
title('Phase Repsonse of Duobinary Singaling')
endfunction
//Result
//-->exec('C:\Users\SENTHILKUMAR\Desktop\Communication_Toolbox\Digital_Communication\New folder\Duobinary_Signaling.sci', -1)
//
//-->[Amplitude_Response,Phase_Response]= Duobinary_Signaling()
//Enter the bit rate= 8Equalizer to compensate aperture effect
function [E]=EqualizerFor_ApertureEff(Ts)
//Equalizer to Compensate for aperture effect
T_Ts = 0.01:0.01:Ts;
E(1) =1;
for i = 2:length(T_Ts)
E(i) = ((%pi/2)*T_Ts(i))/(sin((%pi/2)*T_Ts(i)));
end
a =gca();
a.data_bounds = [0,0.8;0.8,1.2];
plot2d(T_Ts,E,5)
xlabel('Duty cycle T/Ts')
ylabel('1/sinc(0.5(T/Ts))')
title('Normalized equalization (to compensate for aperture effect) plotted versus T/Ts')
endfunctionRepetition code
function [G,H,x]=RepetitionCode(n,k,m)
//Repetition Codes
//n =block of identical 'n' bits
//k =1 one bit
//m = 1;// bit value = 1
I = eye(n-k,n-k);//Identity matrix
P = ones(1,n-k);//coefficient matrix
H = [I P'];//parity-check matrix
G = [P 1];//generator matrix
x = m.*G; //code word
disp(G,'generator matrix');
disp(H,'parity-check matrix');
disp(x,'code word for binary one input');
endfunction
//Result
//-->exec('C:\Users\SENTHILKUMAR\Desktop\Communication_Toolbox\Digital_Communication\RepetitionCode.sci', -1)
//
//-->n = 5
// n =
//
// 5.
//
//-->k =1
// k =
//
// 1.
//
//-->m =1
// m =
//
// 1.
//
//-->[G,H,x]=RepetitionCode(n,k,m)
//
// generator matrix
//
// 1. 1. 1. 1. 1.
//
// parity-check matrix
//
// 1. 0. 0. 0. 1.
// 0. 1. 0. 0. 1.
// 0. 0. 1. 0. 1.
// 0. 0. 0. 1. 1.
//
// code word for binary one input
//
// 1. 1. 1. 1. 1.
// x =
//
// 1. 1. 1. 1. 1.
// H =
//
// 1. 0. 0. 0. 1.
// 0. 1. 0. 0. 1.
// 0. 0. 1. 0. 1.
// 0. 0. 0. 1. 1.
// G =
//
// 1. 1. 1. 1. 1.Reed Solomon Codes
function[n,p,r] = ReedSolomon_Codes(m)
// Reed-Solomon Codes
//Single-error-correcting RS code with a 2-bit byte
//m = m-bit symbol
k = 1^m; //number of message bits
t =1; //single bit error correction
n = 2^m-1; //code word length in 2-bit byte
p = n-k; //parity bits length in 2-bit byte
r = k/n; //code rate
disp(n,'code word length in 2-bit byte n =')
disp(p,'parity bits length in 2-bit byte n-k=')
disp(r,'Code rate:r = k/n =')
disp(2*t,'It can correct any error upto =')
endfunction
//Result
//-->exec('C:\Users\SENTHILKUMAR\Desktop\Communication_Toolbox\Digital_Communication\ReedSolomon_Codes.sci', -1)
//
//-->m=2
// m =
//
// 2.
//
//-->[n,p,r] = ReedSolomon_Codes(m)
//
// code word length in 2-bit byte n =
//
// 3.
//
// parity bits length in 2-bit byte n-k=
//
// 2.
//
// Code rate:r = k/n =
//
// 0.3333333
//
// It can correct any error upto =
//
// 2.
// r =
//
// 0.3333333
// p =
//
// 2.
// n =
//
// 3.
//
