//Caption : FULL Band Echo Cnacellation using NLMS Adaptive Filter
clc;
//Reading a speech signal
[x,Fs,bits]=wavread("E:\4.wav");
order = 40; // Adaptive filter order
x = x';
N = length(x); //length of speech signal
//Delay introduced in echo path
delay = 100;
//Echo at same speaker
xdelayed = zeros(1,N+delay);
for i = delay+1:N+delay
xdelayed(i) = x(i-delay);
end
//Initialize the adaptive filter coefficients to zero
hcap = zeros(1,order);
//To avoid negative values generated during convolution
for i = 1:order-1
xlms(i) = 0.0;
end
for i = order:N+order-1
xlms(i) = x(i-order+1);
end
Power_X = pow_1(x,N); //Average power of speech signal
//Calculation of step size of adaptive filter
delta = 1/(10*order*Power_X);
[x_out,Adapt_Filter_IR] = adapt_filt(xlms,xdelayed,hcap,delta,N,order)
figure(1)
subplot(3,1,1)
plot([1:N],x)
title('Speech Signal Generated by some Speaker A')
sound(x,Fs,16)
subplot(3,1,2)
plot([1:length(xdelayed)],xdelayed)
title('Echo signal of speaker A received by speaker A')
sound(xdelayed,Fs,16)
subplot(3,1,3)
plot([1:length(x_out)],x_out)
title('Echo signal at speaker A after using NLMS Adaptive filter Echo Canceller')
figure(2)
plot2d3('gnn',[1:length(Adapt_Filter_IR)],Adapt_Filter_IR)
title('Adaptive Filter (Echo Canceller) Impulse response')
function [c] = PCM_Encoding(x,L,en_code)
//Encoding: Converting Quantized decimal sample values in to binary
//x = input sequence
//L = number of qunatization levels
//en_code = normalized input sequence
n = log2(L);
c = zeros(length(x),n);
for i = 1:length(x)
for j = n:-1:0
if(fix(en_code(i)/(2^j))==1)
c(i,(n-j)) =1;
en_code(i) = en_code(i)-2^j;
end
end
end
disp(c)
//PROGRAM TO DESIGN AND OBTAIN THE FREQUENCY RESPONSE OF FIR FILTER
//Band PASS FILTER - Window Based
clear all;
clc;
close;
M = 7 //Filter length = 7
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) = (Wc2-Wc1)/%pi;
else
n
hd(n) = (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)
//Caption: Reading a Speech Signal &
//[1]. Write it in another file
//[2]. Changing the bit depth from 16 bits to 8 bits
clear;
clc;
[y,Fs,bits_y]=wavread("E:\4.wav");
wavwrite(y,Fs,8,"E:\4_8bit.wav");
[x,Fs,bits_x]=wavread("E:\4_8bit.wav");
Ny = length(y); //Number of samples in y (4.wav)
Nx = length(x); //Number of samples in x (4_8bit.wav)
Memory_y = Ny*bits_y; //memory requirement for 4.wav in bits
Memory_x = Nx*bits_x; //memory requirement for 4_8bit.wav in bits
disp(Memory_y,'memory requirement for 4.wav in bits =')
disp(Memory_x,'memory requirement for 4_8bit.wav in bits =')
//Result
//memory requirement for 4.wav in bits =
// 133760.
//
//memory requirement for 4_8bit.wav in bits =
// 66880.
function[y,M] = adapt_filt(xlms,B,h,delta,l,l1)
//x = the signal from the speaker directly
//B = the signal thorugh the echo path
//h = impulse response of adaptive filter
//l = length of signal 'x'
//l1 = length of adaptive filter order
for k = 1:150
for n = 1:l
xcap = xlms((n+l1-1):-1:(n+l1-1)-l1+1);
yout(n) = h*xcap;
e(n) = B(n)- yout(n);
xnorm = 0.001+(xcap*xcap');
h = h+((delta*e(n))*(xcap'));
end
eold = 0.0;
for i = 1:l
MSE = eold+(e(i)^2);
eold = MSE;
end
if MSE <= 0.0001 then
break;
end
end
y = zeros(1,length(e));
M = zeros(1,length(h));
y = e;
M = h;
endfunction
//Determination of spectrum of a signal
//With maximum normalized frequency f = 0.4
//using Blackmann window
clear all;
close;
clc;
N = 11;
cfreq = [0.4 0];
[wft,wfm,fr]=wfir('lp',N,cfreq,'re',0);
wft; // Time domain filter coefficients
wfm; // Frequency domain filter values
fr; // Frequency sample points
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));
wft_blmn(n) = wft(n)*h_balckmann(n);
end
wfm_blmn = frmag(wft_blmn,length(fr));
WFM_blmn_dB =20*log10(wfm_blmn);
plot2d(fr,WFM_blmn_dB)
xtitle('Frequency Response of Blackmann window Filtered output N = 11','Frequency in cycles per samples f','Energy density in dB')
//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))
//Using Digital Filter Transformation, the First order
//Analog IIR Butterworth LPF converted into Digital
//Butterworth HPF
clear all;
clc;
close;
s = poly(0,'s');
Omegac = 0.2*%pi;
H = Omegac/(s+Omegac);
T =1;//Sampling period T = 1 Second
z = poly(0,'z');
Hz_LPF = horner(H,(2/T)*((z-1)/(z+1)));
alpha = -(cos((Omegac+Omegac)/2))/(cos((Omegac-Omegac)/2));
HZ_HPF=horner(Hz_LPF,-(z+alpha)/(1+alpha*z))
HW =frmag(HZ_HPF(2),HZ_HPF(3),512);
W = 0:%pi/511:%pi;
plot(W/%pi,HW)
a=gca();
a.thickness = 3;
a.foreground = 1;
a.font_style = 9;
xgrid(1)
xtitle('Magnitude Response of Single pole HPF Filter Cutoff frequency = 0.2*pi','Normalized Digital Frequency W/pi--->','Magnitude');
// 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 -->');
//Design of FIR Filter using Frquency Sampling Technique
//Low Pass Filter Design
//Cutoff Frequency Wc = pi/2
//M = Filter Lenth = 7
clear;
clc;
M = 7;
N = ((M-1)/2)+1
wc = %pi/2;
for k =1:M
w(k) = ((2*%pi)/M)*(k-1);
if (w(k)>=wc)
k-1
break
end
end
for i = 1:k-1
Hr(i) = 1;
G(i) = ((-1)^(i-1))*Hr(i);
end
for i = k:N
Hr(i) = 0;
G(i) = ((-1)^(i-1))*Hr(i);
end
h = zeros(1,M);
for n = 1:M
for k = 2:N
h(n) = G(k)*cos((2*%pi/M)*(k-1)*((n-1)+(1/2)))+h(n);
end
h(n) = (1/M)* (G(1)+2*h(n));
end
[hzm,fr]=frmag(h,256);
hzm_dB = 20*log10(hzm)./max(hzm);
plot(2*fr,hzm_dB)
xtitle('Frequency Response of LPF with Normalized cutoff =0.5','Normalized Frequency W ------>','Magnitude Response H(w)---->');
xgrid(1);
//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)
