//Program for Linear Convolution
clc;
clear all;
close ;
x = input('enter x seq');
h = input('enter h seq');
m = length(x);
n = length(h);
//Method : Using Direct Convolution Sum Formula
for i = 1:n+m-1
conv_sum = 0;
for j = 1:i
if (((i-j+1) <= n)&(j <= m))
conv_sum = conv_sum + x(j)*h(i-j+1);
end;
y(i) = conv_sum;
end;
end;
disp(y,'y=')
// Continuous Time Fourier Transform and Frequency Response of a Square Waveform
// x(t)= A, from -T1 to T1
clear;
clc;
close;
// CTS Signal
A =1; //Amplitude
Dt = 0.005;
T1 = 4; //Time in seconds
t = -T1:Dt:T1;
for i = 1:length(t)
xt(i) = A;
end
//
// Continuous-time Fourier Transform
Wmax = 2*%pi*1; //Analog Frequency = 1Hz
K = 4;
k = 0:(K/1000):K;
W = k*Wmax/K;
xt = xt';
XW = xt* exp(-sqrt(-1)*t'*W) * Dt;
XW_Mag = real(XW);
W = [-mtlb_fliplr(W), W(2:1001)]; // Omega from -Wmax to Wmax
XW_Mag = [mtlb_fliplr(XW_Mag), XW_Mag(2:1001)];
//
subplot(2,1,1);
a = gca();
a.data_bounds=[-4,0;4,2];
a.y_location ="origin";
plot2d2(t,xt);
xlabel('t in msec.');
title('Continuous Time Signal x(t)')
subplot(2,1,2);
a = gca();
a.y_location ="origin";
a.x_location ="origin";
plot(W,XW_Mag);
xlabel('Frequency in Radians/Seconds');
title('Continuous-time Fourier Transform X(jW)')
//Program To convert analog filter into digital filter using
//mapping = (z-(z^-1))/T - Backward Difference
clear all;
clc;
close;
s = poly(0,'s');
H = 1/((s+0.1)^2+9)
T =1;//Sampling period T = 1 Second
z = poly(0,'z');
Hz = horner(H,(1/T)*(z-(z^-1)))
//Program To convert analog IIR filter into digital IIR filter using
//Bilinear Transformation
clear all;
clc;
close;
s = poly(0,'s');
H = (s+0.1)/((s+0.1)^2+16);
T = 0.5;
z = poly(0,'z');
Hz = horner(H,(2/T)*((z-1)/(z+1)))
//Low Pass FIlter of length M = 11
//Pass band Edge frequency fp = 0.1 and a Stop edge frequency fs = 0.15
// Choose the number of cosine functions and create a dense grid
// in [0,0.2) and [0.25,0.5)
//magnitude for pass band = 1 & stop band = 0 (i.e) [1 0]
//Weighting function =[1 1]
clear all;
clc;
close;
M =11;
hn=eqfir(11,[0,0.2;0.25,0.5],[1 0],[1 1]);
[hm,fr]=frmag(hn,256);
disp(hn,'The Filter Coefficients are:')
figure
plot(.5*(0:255)/256,20*log10(frmag(hn,256)));
xlabel('Normalized Digital Frequency fr');
ylabel('Magnitude in dB');
title('Frequency Response of FIR LPF using REMEZ algorithm M=11')
xgrid(2)
//PROGRAM TO DESIGN AND OBTAIN THE FREQUENCY RESPONSE OF FIR FILTER
//HIGH PASS FILTER -WINDOW BASED
clear all;
clc;
close;
M = 7 //Filter length = 7
Wc = %pi/4; //Digital Cutoff frequency
Tuo = (M-1)/2 //Center Value
for n = 1:M
if (n == Tuo+1)
hd(n) = 1-Wc/%pi;
else
hd(n) = (sin(%pi*((n-1)-Tuo)) -sin(Wc*((n-1)-Tuo)))/(((n-1)-Tuo)*%pi);
end
end
//Rectangular Window
for n = 1:M
W(n) = 1;
end
//Windowing Fitler Coefficients
h = hd.*W;
disp('Filter Coefficients are')
h;
[hzm,fr]=frmag(h,256);
hzm_dB = 20*log10(hzm)./max(hzm);
plot(fr,hzm_dB)
xlabel('Normalized Digital Frequency W');
ylabel('Magnitude in dB');
title('Frequency Response 0f FIR HPF using Rectangular window M=7')
//PROGRAM TO DESIGN AND OBTAIN THE FREQUENCY RESPONSE OF FIR FILTER
//Band Stop FILTER (or)Band Reject Filter
clear all;
clc;
close;
M = 7 //Filter length = 11
Wc = [%pi/5,3*%pi/4]; //Digital Cutoff frequency
Wc2 = Wc(2)
Wc1 = Wc(1)
Tuo = (M-1)/2 //Center Value
hd = zeros(1,M);
W = zeros(1,M);
for n = 1:M
if (n == Tuo+1)
hd(n) = 1-((Wc2-Wc1)/%pi);
else hd(n)=(sin(%pi*((n-1)-Tuo))-sin(Wc2*((n-1)-Tuo))+sin(Wc1*((n-1)-Tuo)))/(((n-1)-Tuo)*%pi);
end
if(abs(hd(n))<(0.00001))
hd(n)=0;
end
end
hd
//Rectangular Window
for n = 1:M
W(n) = 1;
end
//Windowing Fitler Coefficients
h = hd.*W;
disp(h,'Filter Coefficients are')
[hzm,fr]=frmag(h,256);
hzm_dB = 20*log10(hzm)./max(hzm);
plot(2*fr,hzm_dB)
xlabel('Normalized Digital Frequency W');
ylabel('Magnitude in dB');
title('Frequency Response 0f FIR BPF using Rectangular window M=7')
xgrid(1)
//Evaluating power spectrum of a discrete sequence
//Using N-point DFT
clear all;
clc;
close;
N =32; //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,64,256,512];
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
title('for frequency separation = 0.06')
subplot(2,2,1)
plot2d3('gnn',k1,abs(Pxx1_fk1))
subplot(2,2,2)
plot2d3('gnn',k2,abs(Pxx1_fk2))
subplot(2,2,3)
plot2d3('gnn',k3,abs(Pxx1_fk3))
subplot(2,2,4)
plot2d3('gnn',k4,abs(Pxx1_fk4))
figure
title('for frequency separation = 0.01')
subplot(2,2,1)
plot2d3('gnn',k1,abs(Pxx2_fk1))
subplot(2,2,2)
plot2d3('gnn',k2,abs(Pxx2_fk2))
subplot(2,2,3)
plot2d3('gnn',k3,abs(Pxx2_fk3))
subplot(2,2,4)
plot2d3('gnn',k4,abs(Pxx2_fk4))
//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');
//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]');
// Continuous Time Fourier Transform and Frequency Response of a Square Waveform
// x(t)= A, from -T1 to T1
clear;
clc;
close;
// CTS Signal
A =1; //Amplitude
Dt = 0.005;
T1 = 4; //Time in seconds
t = -T1:Dt:T1;
for i = 1:length(t)
xt(i) = A;
end
//
// Continuous-time Fourier Transform
Wmax = 2*%pi*1; //Analog Frequency = 1Hz
K = 4;
k = 0:(K/1000):K;
W = k*Wmax/K;
xt = xt';
XW = xt* exp(-sqrt(-1)*t'*W) * Dt;
XW_Mag = real(XW);
W = [-mtlb_fliplr(W), W(2:1001)]; // Omega from -Wmax to Wmax
XW_Mag = [mtlb_fliplr(XW_Mag), XW_Mag(2:1001)];
//
subplot(2,1,1);
a = gca();
a.data_bounds=[-4,0;4,2];
a.y_location ="origin";
plot2d2(t,xt);
xlabel('t in msec.');
title('Continuous Time Signal x(t)')
subplot(2,1,2);
a = gca();
a.y_location ="origin";
a.x_location ="origin";
plot(W,XW_Mag);
xlabel('Frequency in Radians/Seconds');
title('Continuous-time Fourier Transform X(jW)')
//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
//Continuous Time Fourier Series Spectral Coefficients of
//a periodic Cosine signal x(t) = cos(Wot)
clear;
close;
clc;
t = 0:0.01:1;
T = 1;
Wo = 2*%pi/T;
xt = cos(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 0.5049505+1.010D-18i 0 0.5049505-1.010D-18i 0 0 0 0
function [y,X,H] = conv2d2(x,h)
[x1,x2] = size(x);
[h1,h2] = size(h);
X = zeros(x1+h1-1,x2+h2-1);
H = zeros(x1+h1-1,x2+h2-1);
Y = zeros(x1+h1-1,x2+h2-1);
for i = 1:x1
for j = 1:x2
X(i,j)=x(i,j);
end
end
for i =1:h1
for j = 1:h2
H(i,j)=h(i,j);
end
end
disp(X,'x=')
disp(H,'h =')
Y = fft2d(X).*fft2d(H);
disp(Y)
y = ifft2d(Y);
endfunction
/////////////////////////////////////////////
function [a2] = fft2d(a)
//a = any real or complex 2D matrix
//a2 = 2D-DFT of 2D matrix 'a'
m=size(a,1)
n=size(a,2)
// fourier transform along the rows
for i=1:n
a1(:,i)=exp(-2*%i*%pi*(0:m-1)'.*.(0:m-1)/m)*a(:,i)
end
// fourier transform along the columns
for j=1:m
a2temp=exp(-2*%i*%pi*(0:n-1)'.*.(0:n-1)/n)*(a1(j,:)).'
a2(j,:)=a2temp.'
end
for i = 1:m
for j = 1:n
if((abs(real(a2(i,j)))<0.0001)&(abs(imag(a2(i,j)))<0.0001))
a2(i,j)=0;
elseif(abs(real(a2(i,j)))<0.0001)
a2(i,j)= 0+%i*imag(a2(i,j));
elseif(abs(imag(a2(i,j)))<0.0001)
a2(i,j)= real(a2(i,j))+0;
end
end
end
/////////////////////////////////////////////////
function [a] =ifft2d(a2)
//a2 = 2D-DFT of any real or complex 2D matrix
//a = 2D-IDFT of a2
m=size(a2,1)
n=size(a2,2)
//Inverse Fourier transform along the rows
for i=1:n
a1(:,i)=exp(2*%i*%pi*(0:m-1)'.*.(0:m-1)/m)*a2(:,i)
end
//Inverse fourier transform along the columns
for j=1:m
atemp=exp(2*%i*%pi*(0:n-1)'.*.(0:n-1)/n)*(a1(j,:)).'
a(j,:)=atemp.'
end
a = a/(m*n)
a = real(a)
endfunction
//Program for Linear Convolution
clc;
clear all;
close ;
x = input('enter x seq');
h = input('enter h seq');
m = length(x);
n = length(h);
//Method : Using Direct Convolution Sum Formula
for i = 1:n+m-1
conv_sum = 0;
for j = 1:i
if (((i-j+1) <= n)&(j <= m))
conv_sum = conv_sum + x(j)*h(i-j+1);
end;
y(i) = conv_sum;
end;
end;
disp(y,'y=')