Linear Convolution
Direct method using matlab command 'conv'
x = [1 0 0 0 0 0 0 0 2 0 0 0 0 0 0 4];
h = [1 2 3 4 3 2 1];
y = conv(x,h)
script file withought using the command
'conv'
%script file to implement linear convolution
%=============================================
clear all
clc
%input sequences
x1 = [1 0 0 0 0 0 0 0 2 0 0 0 0 0 0 4];
x2 = [1 2 3 4 3 2 1];
if length(x1) > length(x2)
x = x1;h = x2;
else
x = x2;h = x1;
end
% first row in circulant matrix
first_row = h;
first_row = [first_row zeros(1,length(x))];
% generation of circulant matrix
for e = 1:length(first_row)
circMat(e,:) = circshift(first_row,[-1 e-1]);
end
% scaling of row elements
%------------------------------------------------------------------------
for e = 1:length(x)
circMat(e,:) = x(e)*circMat(e,:);
end
circMat = circMat(1:end-length(h),:);
%sum of collumn
%------------------------------------------------------------------------
y = sum(circMat,1);
%display of convoluted signal
%------------------------------------------------------------------------
subplot(3,1,1);stem(x,'linewidth',3);xlim([0 length(x) + length(h) + 1]);title('input
sequence');grid on
subplot(3,1,2);stem(h,'linewidth',3);xlim([0 length(x) + length(h) + 1]);title('system
function');grid on
subplot(3,1,3);stem(y,'linewidth',3);xlim([0 length(x) + length(h) + 1]);title('system
output');grid on
Linear Convolution using DFT
%script file to implement linear convolution using DFT
%==============================================================
============
clear all
clc
%input sequences
%--------------------------------------------------------------------------
x = [1 0 0 0 0 0 0 0 2 0 0 0 0 0 0 4];
h = [1 2 3 4 3 2 1];
%deciding the length of sequence for DFT
%--------------------------------------------------------------------------
tot_length = length(x) + length(h) + 1;
n = ceil(log(tot_length)/log(2));
N = 2^n;
%zero padding to make the sequence length equal
%--------------------------------------------------------------------------
x(end+1:N) = 0;
h(end+1:N) = 0;
%fft
X = fft(x,N);
H = fft(h,N);
Y = X.*H;
y = ifft(Y);
%display of signals
%--------------------------------------------------------------------------
subplot(3,1,1);stem(x,'linewidth',3);xlim([0 length(x) + length(h) + 1]);title('input
sequence');grid on
subplot(3,1,2);stem(h,'linewidth',3);xlim([0 length(x) + length(h) + 1]);title('system
function');grid on
subplot(3,1,3);stem(y,'linewidth',3);xlim([0 length(x) + length(h) + 1]);title('system
output');grid on
Circular Convolution
%script file to implement circular convolution
%==============================================================
============
clear all
clc
%input sequences
%--------------------------------------------------------------------------
x = [1 0 0 0 0 0 0 0 2 0 0 0 0 0 0 4];
h = [1 2 3 4 3 2 1];
%zero padding to get signals of equal length
%--------------------------------------------------------------------------
if length(x) > length(h)
h(end+1 : length(x)) = 0;
elseif length(h) > length(x)
h(end+1 : length(h)) = 0;
end
%circulant matrix generation
for e = 1:length(h)
circMat(e,:) = circshift(h,[-1 e-1]);
end
for e = 1:length(x)
circMat(e,:) = x(e)*circMat(e,:);
end
%obtaining the o/p signal
%--------------------------------------------------------------------------
y = sum(circMat,1);
%display of convoluted signal
%--------------------------------------------------------------------------
subplot(3,1,1);stem(x,'linewidth',3);title('input sequence');grid on
subplot(3,1,2);stem(h,'linewidth',3);title('system function');grid on
subplot(3,1,3);stem(y,'linewidth',3);title('circular convoluted signal');grid on
Circular Convolution using DFT
%script file to implement circular convolution using DFT
%==============================================================
============
clear all
clc
%input sequences
%--------------------------------------------------------------------------
x = [1 0 0 0 0 0 0 0 2 0 0 0 0 0 0 4];
h = [1 2 3 4 3 2 1];
%deciding the length of sequence for DFT
%--------------------------------------------------------------------------
tot_length = max(length(x),length(h));
n = ceil(log(tot_length)/log(2));
N = 2^n;
%zero padding to make the sequence length equal
%--------------------------------------------------------------------------
x(end+1:N) = 0;
h(end+1:N) = 0;
%fft
%--------------------------------------------------------------------------
X = fft(x,N);
H = fft(h,N);
Y = X.*H;
y = ifft(Y);
%display of convoluted signal
%--------------------------------------------------------------------------
subplot(3,1,1);stem(x,'linewidth',3);title('input sequence');grid on
subplot(3,1,2);stem(h,'linewidth',3);title('system function');grid on
subplot(3,1,3);stem(y,'linewidth',3);title('circular convoluted signal');grid on
Sampling Theorem
%% sampling theorem verification
%==============================================================
============
clear all
close all
clc
%% defining sampled signal
% Given continuous signal has the frequency 1Hz
% sampling rate 'samp_rate' is specified as a multiple of maximum frequency
% if samp_rate >= 2, sampling takes place above nyquist rate and signal
% can be reconstructed from the samples
% if samp_rate < 2, sampling takes place below nyquist rate and signal
% cant be reconstructed becouse of aliasing
samp_rate = 5;
%==============================================================
============
t = 0:(1/300):20; %time definition, sampling frequency of the filter =
300Hz
sa = cos(2*pi*1*t); %approximate analog signal
for e = 1:length(t)
if rem(e,ceil(300/samp_rate)) == 0 %sampling
s(e) = cos(2*pi*1*t(e));
else
s(e) = 0;
end
end
%% disign of lowpass filter(using filter design toolbox)
%==============================================================
============
Fs = 300; % Sampling Frequency
Fpass = 1; % Passband Frequency
Fstop = 3; % Stopband Frequency
Dpass = 0.057501127785; % Passband Ripple
Dstop = 0.0001; % Stopband Attenuation
dens = 20; % Density Factor
% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fpass, Fstop]/(Fs/2), [1 0], [Dpass, Dstop]);
% Calculate the coefficients using the FIRPM function.
b = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);
%% filtering the descrete sequence
%==============================================================
============
num_coeff = Hd.Numerator;
den_coeff = [1];
initial_cond = zeros(1,max(length(num_coeff),length(den_coeff)) - 1);
lp_s = filter(num_coeff,den_coeff,s,initial_cond);
%% display of the signals
subplot(2,1,1);
plot(t(round(end/2):end),sa(round(end/2):end),'g','linewidth',2);;hold on
stem(t(round(end/2):end),s(round(end/2):end),'linewidth',2);
title('sampled signal');
subplot(2,1,2);
plot(t(round(end/2):end),lp_s(round(end/2):end),'r','linewidth',2);
title('reconstructed signal');
Solution of Difference Equation
ex - for the system described by defference equation
y[n] = x[n] - x[n-1] ------(1)
compare the system with standard equation
a1*y[n] + a2*y[n-1] + a3*y[n-2]... = b1*x[n] + b2*x[n-1] + b3*x[n-2]... ------(2)
a = [a1 a2 a3 ....] = [1 0 0 ...]
b = [b1 b2 b3 ....] = [1 -1 0 ...]
%% program to solve difference equation
%% input
x = [1 2 3 4 5 6 7 8 9 10]
%% previous input and output
xi = [0];
yi = [0];
%% coefficients of difference equation
a = [1 0 0 0];
b = [1 -1 0 0];
%% finding the initial condition
initial_condition = filtic(b,a,xi,yi);
%% forced response of the system (yf)
%o/p because of i/p under zero initial condition
yf = filter(b,a,x,zeros(size(initial_condition)));
%% natural response
%o/p because of initial condition when input is zero
yn = filter(b,a,zeros(size(x)),initial_condition);
%% total response
yt = yn + yf;
%% display of result
subplot(2,2,1); stem(x,'linewidth',3); title('input'); grid on
subplot(2,2,2); stem(yt,'linewidth',3); title('total response'); grid on
subplot(2,2,3); stem(yf,'r','linewidth',3); title('forced response'); grid on
subplot(2,2,4); stem(yn,'r','linewidth',3); title('natural response'); grid on
Impulse Response of a given System
Impulse response is output of the system when input is unit impulse
ex - for the system described by defference equation
y[n] + 0.6y[n-1] + 0.6y[n-2]= x[n] ------(1)
compare the system with standard equation
a1*y[n] + a2*y[n-1] + a3*y[n-2]... = b1*x[n] + b2*x[n-1] + b3*x[n-2]... ------(2)
a = [a1 a2 a3 ....] = [1 0.6 0.6 0 0 0 ...]
b = [b1 b2 b3 ....] = [1 0 0 0 0 0...]
take input as x = [1 0 0 0 0 0 0 ...]
%% program to find impuse response
%% input
x = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
%% coefficients of difference equation
a = [1 0.6 0.6 0];
b = [1 0 0];
%% Impulse response of the system (yf)
%o/p because of i/p under zero initial condition
Imp_resp = filter(b,a,x);
%% display of result
stem(Imp_resp,'linewidth',3); title('impulse response '); grid on
Auto correlation
correlation of signal with itself is called Auto correlation.
1. always auto correlation has maximum at 0. this is called Cauchy–Schwarz inequality.
2. Rxx(n) < Rxx(0) i.e all the values of Rxx are either less than ar equal to Rxx(0)[value
at origin]
3. maximum value of auto correlation function is equal to energy of the given signal
4. DFT(Rxx) = DFT(x)*conjugate(DFT(x))
clear all
close all
clc
%==============================================================
%% i/p sequence
x = [1 2 3 4];
Rxx = conv(x,fliplr(x)); %correlation operation
x_length = length(x); %x_length = length of the sequence x
N = length(Rxx); %length of Rxx
n = -(x_length-1):(x_length-1); %finding the values of n for display of the signal
%==============================================================
%% property-1 maximum of auto correlation is at origin
disp('property-1');
[max_value max_index] = max(Rxx); %find the maximum value and its index
if (max_value == Rxx(x_length))
disp('maximum is at origin');
end
%==============================================================
%% property-2 Rxx(n) <= Rxx(0)
disp('property-2');
if sum(Rxx > Rxx(x_length)) == 0 %check for the element in Rxx > max_value i.e
Rxx at origin
disp('Rxx(n) <= Rxx(0)');
end
%==============================================================
%% property-3 value of Rxx at origin = total energy of the given signal
disp('property-3');
if Rxx(x_length) == sum(x.^2)
disp('value of Rxx at origin = total energy of the given signal');
end
%==============================================================
%% property-4 DFT(Rxx) = DFT(x)*conjugate(DFT(x))
disp('property-4');
LHS = abs(fft(Rxx,N));
RHS = fft(x,N).*conj(fft(x,N));
LHS = uint8(LHS);
RHS = uint8(RHS);
if LHS == RHS
disp('DFT(Rxx) = DFT(x)*conjugate(DFT(x))');
end
%==============================================================
%% plot of signals
figure, %i/p sequence
stem(n,[zeros(1,x_length-1),x]);
title('Input sequence');
figure, %prop-1,prop-2 and prop-3
subplot(2,1,1);
stem(n,Rxx);
title('Auto-correlation function');
subplot(2,1,2);
stem(n,[zeros(1,x_length-1),x.^2]);
title('energy of sequence x');
figure, %prop-4
subplot(2,1,1);
stem(LHS);
title('DFT of Rxx');
subplot(2,1,2);
stem(RHS);
title('(DFT of x)*conj(DFT of x)');
Cross-correlation
clear all
close all
clc
%% i/p sequence
x = [1 2 3 4 5 6];
y = [1 2 3]
Rxy = conv(x,fliplr(y)); %cross-correlation operation by time reversal of y
Ryx = conv(y,fliplr(x)); %cross-correlation operation by time reversal of x
x_length = length(x); %x_length = length of the sequence x
y_length = length(y); %y_length = length of the sequence y
N = length(Rxy); %length of Rxy = length of Ryx
%==============================================================
=======
%% property-1 cross correlation is not commutative
disp('property-1');
if Rxy ~= Ryx
disp('cross correlation is not commutative');
end
%==============================================================
=======
%% property-2 Rxy(n) = Ryx(-n)
disp('property-2');
if Rxy == fliplr(Ryx)
disp('Rxy(n) = Ryx(-n)');
end
%==============================================================
=======
%% property-3 DFT(Rxy) = DFT(x)*conjugate(DFT(y))
disp('property-3');
LHS = abs(fft(Rxy,N));
RHS = abs(fft(x,N).*conj(fft(y,N)));
LHS = uint8(LHS);
RHS = uint8(RHS);
if LHS == RHS
disp('DFT(Rxx) = DFT(x)*conjugate(DFT(y))');
end
%==============================================================
=======
%% plot of signals
subplot(2,3,1); %i/p signal-x
stem(x);
title('sequence - x');
subplot(2,3,4); %i/p signal-y
stem(y);
title('sequence - y');
subplot(2,3,2); %Rxy
stem(Rxy);
title('Rxy');
subplot(2,3,5); %Ryx
stem(Ryx);
title('Ryx');
subplot(2,3,3); %DFT(Rxx)
stem(LHS);
title('DFT(Rxx)');
subplot(2,3,6); %DFT(x)*conjugate(DFT(x))
stem(RHS);
title('DFT(x)*conjugate(DFT(x))');
Circular Convolution (tested OK)
#include<stdio.h>
void main()
{
float x[100],h[100],sum= 0.0,y[100];
int m,n,len1,len2,i,k;
printf("enter the length of the sequence x=");
scanf("%d",&len1);
printf("Enter the length of the sequence h=");
scanf("%d",&len2);
if(len1==len2)
{
printf("\n Enter the values of the sequnce x=");
for(i=0;i<len1;i++)
{
printf("\n x[%d]=",i);
scanf("%f",&x[i]);
printf("\t%f",x[i]);
}
printf("\nEnter the value of sequence h=");
for(i=0;i<len1;i++)
{
printf("\nh[%d]=",i);
scanf("%f",&h[i]);
printf("\t %f",h[i]);
}
printf("\n After convolution =");
for(n=0;n<=len1;n++)
{
sum=0.0;
for(m=0;m<len2;m++)
{
if((n-m)<0)
k=n-m+len2;
else
k=(n-m);
sum=sum+(x[m]*h[k]);
}
y[n]=sum;
printf("\ny[%d]=%f",n,y[n]);
}
}
else
printf("\ncircular convolution not possible\n");
}
sample o/p
enter the length of the sequence x=4
Enter the length of the sequence h=4
Enter the values of the sequnce x=
x[0]=1
1.000000
x[1]=2
2.000000
x[2]=3
3.000000
x[3]=4
4.000000
Enter the value of sequence h=
h[0]=1
1.000000
h[1]=2
2.000000
h[2]=3
3.000000
h[3]=4
4.000000
After convolution =
y[0]=26.000000
y[1]=28.000000
y[2]=26.000000
y[3]=20.000000
Impulse Response (tested OK)
#include<stdio.h>
void main()
{
int i,n,L;
float a[3],b[3],x[100],y[100];
printf("enter the length of the output sequence=");
scanf("%d",&L);
printf("\n enter the numerator coefficients(or) coefficients of'x'=");
for (i=0;i<3;i++)
{
printf("\nb[%d]=",i);
scanf("%f",&b[i]);
}
printf("\n enter the denominator coefficients(or) coeff of 'y'=");
for (i=0;i<3;i++)
{
printf("\na[%d]=",i);
scanf("%f",&a[i]);
}
x[0]=1;
for (n=1;i<L;i++)
{
x[n]=0;
y[0]=b[0]*x[0]/a[0];
y[1]=((b[1]*x[0])-(a[1]*y[0]))/a[0];
}
for(n=2;n<=L;n++)
{
y[n]=((b[2]*x[n-2])-(a[1]*y[n-1])-(a[2]*y[n-2]))/a[0];
}
printf("\n impulse response=");
for (i=0;i<L;i++)
{
printf("\n y[%d]=%f",i,y[i]);
}
}
sample output
enter the length of the output sequence=10
enter the numerator coefficients(or) coefficients of'x'=
b[0]=1
b[1]=2
b[2]=3
enter the denominator coefficients(or) coeff of 'y'=
a[0]=1
a[1]=-0.2
a[2]=0
impulse response=
y[0]=1.000000
y[1]=2.200000
y[2]=3.440000
y[3]=0.688000
y[4]=000000000000000000000000.000000
y[5]=0.000000
y[6]=000000000000000000000000.000000
y[7]=0.000000
y[8]=000000000000000000000000.000000
DFT (tested OK)
#include<stdio.h>
#include<math.h>
void main()
{
int n,k,i,l;
float x[100],x_real[100],x_img[100],x_phase[100],sum_r=0.0,sum_i=0.0,pi,x_mag[100];
printf("\n Enter the length of sequence =");
scanf("%d",&l);
printf("\n Enter the element of the sequence x");
for(i=0;i<l;i++)
{
printf("\n x[%d]=",i);
scanf("%f",&x[i]);
}
pi=3.14159265358979;
for(k=0;k<l;k++)
{
for(n=0;n<l;n++)
{
sum_r=sum_r+( x[n]*cos(2*pi*n*k/l));
sum_i=sum_i+( x[n]*sin(2*pi*n*k/l));
}
x_real[k]=sum_r;
x_img[k]=sum_i;
sum_i=0.0;
sum_r=0,0;
}
printf("\n DFT of the sequence=\n");
for(i=0;i<l;i++)
printf("X[%d]=%f+i%f\n",i,x_real[i],x_img[i]);
printf("\n magnitude of DFT=\n");
for(i=0;i<n;i++)
{
x_mag[i]=sqrt((x_real[i]*x_real[i])+(x_img[i]*x_img[i]));
printf("\n |X[%d]|=%f\n",i,x_mag[i]);
}
printf("\n phase=\n");
for(i=0;i<l;i++)
{
x_phase[i]=x_img[i]/x_real[i];
printf("\n Phase (X[%d])=arctan(%f)\n",i,x_phase[i]);
}
}
Sample output
Enter the length of sequence =4
Enter the element of the sequence x
x[0]=1
x[1]=2
x[2]=3
x[3]=4
DFT of the sequence=
X[0]=10.000000+i0.000000
X[1]=-2.000000+i-2.000000
X[2]=-2.000000+i0.000000
X[3]=-2.000000+i2.000000
magnitude of DFT=
|X[0]|=9.999999
|X[1]|=2.828428
|X[2]|=2.000000
|X[3]|=2.828427
