app4/Notch_Filter.m at master · geappteam/app4 · GitHub

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
clc
close all
clear all %#ok<CLSCR>

% Problematique
[y,Fs] = audioread('res/note_basson_plus_sinus_1000_Hz.wav');

% ***
% Bandpass Filter : 2 pole 2 zeros
% Dc Gain = 0dB
% Cut-off frequency between 980 - 1020 Hz
%                        -> Must be over -3 dB
% ***

% For information : "Check digitalFilter" class with Tool Help !

%% Bandpass Filter (IIR)
r = 0.9985;
N = length(y);
t = (0:N-1)'/Fs;                                    % t : time axis
fn_Low = 980;                                       % 980 Hz
fn_High = 1020;                                     % 1020 Hz
fn = ((fn_Low + fn_High)/2);                        % 1000 Hz

%% Put poles at same angle, inside unit circle
b = [1 -2*cos(2*pi*fn/Fs) 1];                       % Filter Coefficients
a = [1 -2*r*cos(2*pi*fn/Fs) r^2];                   % Filter Coefficients
b = b/sum(b)*sum(a);                                % Make Dc gain to 1

[z,p] = tf2zp(b,a);                                 % Check pole and zero

[H,w] = freqz(b,a);                                 % Freqz returns the frequency response based
                                                    % on the current filter coefficients
y_Out = filter(b, a, y);                            % Apply BandpassFilter

% Check Transfer Function
sys = filt(b,a)                                     % Form (z^-1)
sys_Z = tf(b,a,0.1)                                 % Form (z) Ts = 0.1

% Converting amplitude to dB amplitude
H_dB = 20*log10(H);

%% Display IIR pole-zero diagram
figure('Name','Display IIR pole-zero diagram')
clf
zplane(b,a)
title('IIR notch filter')

%% Display IIR Frequency response
figure('Name','Display IIR Frequency response')
clf
plot(w/(2*pi)*Fs, H_dB)
xlabel('Frequency (Hz)')
title('Frequency response (IIR) filter')
box off

%% Display IIR filter to noisy speech
figure('Name','Display IIR filter to noisy speech')
clf
plot(t,y,t,y_Out)
legend('Noisy signal','Filtered (IIR)')
xlabel('Time (sec)')

%% Display IIR filter impulse response
figure('Name','IIR impulse response')
clf
impz(b,a,44100)

%% Reponse a un signal de 1khz
figure('Name','IIR: Reponse sinus 1khz')
clf
td = 0:8000;
tt = td / Fs;
signal  = sin(2*pi*1000*tt);
plot(tt, filter (b,a,signal));
ylabel('Amplitude')
xlabel('Temps (secondes)')
title('IIR: reponse sin (1khz)')