-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfir_filter.m
152 lines (129 loc) · 4.02 KB
/
fir_filter.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
%FIR Bandstop Filter
close all;
clear all;
indexNo = 150009;
%Get A, B, C values from the index number
A = fix(mod(indexNo,1000)/100);
B = fix(mod(indexNo,100)/10);
C = mod(indexNo,10);
%Define filter specifications
Ap = 0.05 + 0.01*A; %Passband ripple
Aa = 40 + B; %Minimum stopband attenuation
Omega_p1 = C*100 + 300; %Lower passband frequency
Omega_p2 = C*100 + 850; %Upper passband frequency
Omega_a1 = C*100 + 400; %Lower stopband frequency
Omega_a2 = C*100 + 700; %Upper stopband frequency
Omega_s = 2*(C*100 + 1200); %Sampling frequency
T = 2*pi/Omega_s;
Bt = min((Omega_a1 - Omega_p1),(Omega_p2 - Omega_a2)); %Transition width
Omega_c1 = Omega_p1 + Bt*0.5; %Lower cut-off frequency
Omega_c2 = Omega_p2 - Bt*0.5; %Upper cut-off frequency
%Choose delta value
delta_p = (10^(0.05*Ap)-1)/(10^(0.05*Ap)+1);
delta_a = 10^(-0.05*Aa);
delta = min(delta_p, delta_a);
%Actual Stopband attenuation
Aaa = -20*log10(delta);
%Choose parameter alpha
if Aaa<=21
alpha = 0;
elseif (21<=Aaa) && (Aaa<=50)
alpha = 0.5842*(Aaa-21)^0.4 + 0.07886*(Aaa-21);
elseif Aaa>50
alpha = 0.1102*(Aaa-8.7);
end
%Choose parameter D
if Aaa<=21
D = 0.9222;
else
D = (Aaa-7.95)/14.36;
end
%Choose the lowest odd value of N
if mod(ceil((Omega_s*D)/Bt+1),2) == 0;
N = ceil((Omega_s*D)/Bt+1) + 1;
else
N = ceil((Omega_s*D)/Bt+1);
end
%Plot the Window function (wk) from Kaiser window function
nr = -(N-1)/2 : 1 : (N-1)/2; %define the range where wk is non-zero.
beta = alpha*(1 - ((2*nr)/(N-1)).^2).^0.5;
I_beta = 1; I_alpha = 1;
for k = 1 : 1 : 100
I_beta = I_beta + ((1/factorial(k))*(beta/2).^k).^2;
I_alpha = I_alpha + ((1/factorial(k))*(alpha/2).^k).^2;
end
wk = I_beta./I_alpha;
figure;
stem(nr,wk);
title('Windowing function');
xlabel('n'); ylabel('w[n]');
grid on;
%Compute h[n]
n1 = -(N-1)/2 : 1 : -1; %Range for negative values
n2 = 1 : 1 : (N-1)/2; %Range for positive values
h1 = (((1/pi)./n1).*(sin(Omega_c1*n1*T) - sin(Omega_c2*n1*T)));
h2 = (((1/pi)./n2).*(sin(Omega_c1*n2*T) - sin(Omega_c2*n2*T)));
h0 = 1 + (2*(Omega_c1-Omega_c2))/Omega_s;
n = [n1, 0, n2];
hn = [h1, h0, h2]; %h[n] array
figure;
stem(n,hn);
title('Impulse Response of Ideal Bandstop Filter');
xlabel('n'); ylabel('h[n]');
grid on;
%Compute Digital filter
filt = hn.*wk;
figure;
stem(n, filt);
title('Impulse Response of Non-causal Filter');
xlabel('n'); ylabel('H');
grid on;
figure;
n_new = 0:1:(N-1);
stem(n_new, filt);
title('Impulse Response of Causal Filter');
xlabel('n'); ylabel('H');
grid on;
%Magnitude Response
fvtool(filt)
%compute the Omega values and plot the Excitation in time domain
w1 = Omega_p1/2;
w2 = (Omega_a1+Omega_a2)/2;
w3 = (Omega_p2+Omega_s/2)/2;
ns = 0:1:300; %No. of samples
xnT = sin(w1*ns*T)+sin(w2*ns*T)+sin(w3*ns*T); %Excitation function
figure;
plot(ns, xnT);
title('Excitation');
xlabel('time(s)');
ylabel('x(nT)');
y = conv2(xnT,filt);
figure;
plot([1:length(y)]*T*(length(xnT))/(length(y)),y);
title('Filtered Signal');
xlabel('time(s)');
N_FFT = 2^nextpow2(numel(ns)); %Next power of 2 from length of y
xnT_FFT = fft (xnT, N_FFT)/numel(ns);
f = (Omega_s)/2* linspace(0, 1, N_FFT/2+1);
figure;
plot(f, 2*abs(xnT_FFT(1:N_FFT/2+1)));
title('Discrete Fourier Transform Excitation');
xlabel('Frequency(rad/s)');
N_FFT = 2^nextpow2(numel(ns)); %Next power of 2 from length of y
Y_FFT = fft(y, N_FFT)/numel(ns);
f = (Omega_s)/2*linspace(0, 1, N_FFT/2+1);
figure;
plot(f, 2*abs(Y_FFT(1:N_FFT/2+1)));
title('Discrete Fourier Transform Filtered from Kaiser Window Bandstop Filter');
xlabel('Frequency(rad/s)');
%Ideal Filter
xnT=sin(w1*ns*T)+sin(w3*ns*T);
Y=conv2(xnT,filt);
N_FFT = 2^nextpow2(numel(ns)); %Next power of 2 from length of y
Y_FFT = fft(Y, N_FFT)/numel(ns);
f = (Omega_s)/2*linspace(0, 1, N_FFT/2+1);
figure;
plot(f, 2*abs(Y_FFT(1:N_FFT/2+1)));
title('DFT Filtered from Ideal Bandstop Filter');
xlabel('Frequency(rad/s)');
%%%%%%%%%%%%%end%%%%%%%%%%%%%