## Spectrum of signal (Frequency Response)- Blackmann Window

March 24, 2011 Coded in Scilab
``````//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')``````

## Power Spectrum using N-point DFT

March 24, 2011 Coded in Scilab
``````//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))``````

## Digital Filter Transformation - LPF to BPF (correct)

March 24, 2011 Coded in Scilab
``````//Transformation of  LPF IIR Digital Butterworth Filter into BPF
clear all;
clc;
close;
omegaP = 0.2*%pi;
omegaL =  (2/5)*%pi;
omegaU =  (3/5)*%pi;
z=poly(0,'z');
H_LPF = (0.245)*(1+(z^-1))/(1-0.509*(z^-1))
alpha = (cos((omegaU+omegaL)/2)/cos((omegaU-omegaL)/2));
k = (cos((omegaU - omegaL)/2)/sin((omegaU - omegaL)/2))*tan(omegaP/2);
NUM =-((z^2)-((2*alpha*k/(k+1))*z)+((k-1)/(k+1)));
DEN = (1-((2*alpha*k/(k+1))*z)+(((k-1)/(k+1))*(z^2)));
HZ_BPF=horner(H_LPF,NUM/DEN)
disp(HZ_BPF,'Digital BPF IIR Filter H(Z)= ')
HW  =frmag(HZ_BPF(2),HZ_BPF(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 BPF Filter cutoff frequency [0.4,0.6]','Normalized Digital Frequency--->','Magnitude');``````

## Digital Filter Transformation - LPF to BSF

March 24, 2011 Coded in Scilab
``````//Transformation of  LPF IIR Digital Butterworth Filter into BSF
clear all;
clc;
close;
omegaP = 0.2*%pi;
omegaL =  (2/5)*%pi;
omegaU =  (3/5)*%pi;
z=poly(0,'z');
H_LPF = (0.245)*(1+(z^-1))/(1-0.509*(z^-1))
alpha = (cos((omegaU+omegaL)/2)/cos((omegaU-omegaL)/2));
k = tan((omegaU - omegaL)/2)*tan(omegaP/2);
NUM =((z^2)-((2*alpha/(1+k))*z)+((1-k)/(1+k)));
DEN = (1-((2*alpha/(1+k))*z)+(((1-k)/(1+k))*(z^2)));
HZ_BPF=horner(H_LPF,NUM/DEN)
HW  =frmag(HZ_BPF(2),HZ_BPF(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 BSF Filter cutoff freq [0.4,0.6] ','Normalized Digital Frequency--->','Magnitude');``````

## Digital Filter Transformation - LPF to BPF

March 24, 2011 Coded in Scilab
``````//Transformation of  LPF IIR Digital Butterworth Filter into BSF
clear all;
clc;
close;
omegaP = 0.2*%pi;
omegaL =  (2/5)*%pi;
omegaU =  (3/5)*%pi;
z=poly(0,'z');
H_LPF = (0.245)*(1+(z^-1))/(1-0.509*(z^-1))
alpha = (cos((omegaU+omegaL)/2)/cos((omegaU-omegaL)/2));
k = tan((omegaU - omegaL)/2)*tan(omegaP/2);
NUM =((z^2)-((2*alpha/(1+k))*z)+((1-k)/(1+k)));
DEN = (1-((2*alpha/(1+k))*z)+(((1-k)/(1+k))*(z^2)));
HZ_BPF=horner(H_LPF,NUM/DEN)
HW  =frmag(HZ_BPF(2),HZ_BPF(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 BSF Filter cutoff freq [0.4,0.6] ','Normalized Digital Frequency--->','Magnitude');``````

## Digital IIR Filter Transformatin - LPF to HPF

March 24, 2011 Coded in Scilab
``````//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');``````

## Analog Filter transformations -IIR Butterworth Filter

March 24, 2011 Coded in Scilab
``````// To Convert Analog LPF into .High Pass .Band Pass IIR Butterworth Filter
//Using Analog Filter Transformations
//For the given cutoff frequency Wc = 500 Hz
clear all;
clc;
close;
omegap =  500;
omegas =  1000;
delta1_in_dB = -3;
delta2_in_dB = -40;
delta1 = 10^(delta1_in_dB/20)
delta2 = 10^(delta2_in_dB/20)
//Calculation of Filter Order
N = log10((1/(delta2^2))-1)/(2*log10(omegas/omegap))
N = ceil(N)
omegac = omegap;
//Poles and Gain Calculation
[pols,gain]=zpbutt(N,omegac);
N =1;
//
omega_LPF = omegap;  //Analog LPF Cutoff frequency
omega_HPF = omega_LPF;  //Analog HPF Cutoff frequency
omega2 = 600;    //Upper Cutoff frequency
omega1 = 300;     //Lower Cutoff Frequency
omega0 = (omega2*omega1);
BW =  omega2 -  omega1;  //Bandwidth
disp('Analog LPF Transfer function')
[hs,pols,zers,gain] = analpf(N,'butt',[0,0],omega_LPF)
hs_LPF = hs;
hs_LPF(2) = hs_LPF(2)/500;
hs_LPF(3)= hs_LPF(3)/500;
s =poly(0,'s');
disp('Analog HPF Transfer function')
h_HPF = horner(hs_LPF,omega_LPF*omega_HPF/s)
disp('Analog BPF Transfer function')
num = (s^2)+omega0
den = BW*s
h_BPF = horner(hs_LPF,omega_LPF*(num/den))
//Plotting Low Pass Filter Frequency Response
figure
fr=0:.1:1000;
hf=freq(hs_LPF(2),hs_LPF(3),%i*fr);
hm=abs(hf);
plot(fr,hm)
a=gca();
a.thickness = 3;
a.foreground = 1;
a.font_style = 9;
xgrid(1)
xtitle('Magnitude Response of LPF Filter Cutoff frequency = 500Hz','Analog Frequency--->','Magnitude');
//Plotting High Pass Filter Frequency Response
figure
fr=0:.1:1000;
hf_HPF=freq(h_HPF(2),h_HPF(3),%i*fr);
hm_HPF=abs(hf_HPF);
plot(fr,hm_HPF)
a=gca();
a.thickness = 3;
a.foreground = 1;
a.font_style = 9;
xgrid(1)
xtitle('Magnitude Response of HPF Filter Cutoff frequency = 500Hz','Analog Frequency--->','Magnitude');
//Plotting Band Pass Filter Frequency Response
figure
fr=0:.1:1000;
hf_BPF=freq(h_BPF(2),h_BPF(3),%i*fr);
hm_BPF=abs(hf_BPF);
plot(fr,hm_BPF)
a=gca();
a.thickness = 3;
a.foreground = 1;
a.font_style = 9;
xgrid(1)
xtitle('Magnitude Response of BPF Filter Upper Cutoff frequency = 600Hz & Lower Cutoff frequency = 300Hz','Analog  Frequency--->','Magnitude');``````

## Analog IIR Butterworth Filter - Scilab

March 22, 20111 comment Coded in 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 -->');``````

## Design of FIR Filter using Frquency Sampling Technique-LPF

March 22, 2011 Coded in Scilab
``````//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);``````

## FIR - BSF - Window Based

March 22, 2011 Coded in Scilab
``````//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)``````