Passing audio file through a notch filter transfer function

I am looking to pass an audio signal through a notch filter to attenuate 1000Hz and 5000Hz frequencies. Here is my code:

[SigN, Fs] = audioread('\\ucofs\NoBackups\UCO VDI\Student Folder Redirection\nblair1\Desktop\Matlab1\RainFireNoise.wav');

Xs=fft(SigN);

L = length(Xs);

t = 0:1:length(Xs)-1;

k = 1;

f = 0:0.1:6000;

w = 2*pi*f;

wo = 2*pi*1000;

w2 = 2*pi*5000;

theta = 85*pi/180;

% zeros and poles for notch filter of frequencies 1000hz and 5000hz

z = [wo*1i; -wo*1i; w2*1i; -w2*1i];

p = [-wo*cos(theta)+wo*sin(theta)*1i; -wo*cos(theta)-wo*sin(theta)*1i; -w2*cos(theta)-w2*sin(theta)*1i; -w2*cos(theta)+w2*sin(theta)*1i];

[NumCof, DenCof] = zp2tf(z,p,k);

Hs = tf(NumCof, DenCof);

Ys = lsim(Hs, Xs, t);

plot(Fs/L*(0:L-1), abs(Ys));

xlim([0 6000]);

My thinking is that I can convert my audio file (SigN) to the laplace domain using the fft (Xs), and then pass it through my notch filter (Hs) using lsim. I know there are built in functions to do this but I'd like to do it using the input signal and transfer function.

The problem I am having is that Hs isnt affecting the frequency components at all and I am unsure why. I am also new to learning about this and am clueless. Any help would be appreciated, thanks.

Answer 1

[William Rose]

@Nathaniel,

[edit: remove junk text at the end]

lsim() wants the input to be the time domain signal. You are passing the frequency domain signal Xs to lsim(). Maybe fixing that will help. I encourage you to think about the frequencies you are using to define the pole and zero locations, and the time vector you pass to lsim(). You define the pole-zero locations with frequencies of 1000 and 5000 Hz. But the time vector, t, which you pass to lsim() is [0,1,2,3,...], which implies a sampling rate of 1 Hz. A 1 Hz-sampled signal does not cotain power at any frequencies above 0.5 Hz. To address the mismatch between the implied sampling rate and the desired filter frequencies, adjust the time vector, which should be t=[0,1,2,...]/fs; where fs=sampling rate in Hz. You define vectors f and w, but you do not do anything with them. What was your goal for those vectors?

Your code and your written description of hte notch filter mention two frequencies, 1 and 5 kHz. Do you want two sharp notch filters - one at 1 kHz and one at 5 kHz? Or do you want a "wide" notch (i.e. a band-stop) covering 1 to 5 kHz? Your poles and zeros suggest two sharp notches, but you only mention one notch filter in your text.

I am not sure why you want to use lsim(). I wouldn't. I would start with a standard filter routine. If it produces unacceptable resuts - and that can happen, especially with higher order filters - then try something else. Here is a notch filter at 1 kHz, designed and then applied to the signal

x=load('mysignal.txt');     % load audio signal file saved as text file
fs=44.1e3;                  % specify sampling rate (Hz)
bsFilter = designfilt('bandstopfir', ...       % response type
    'FilterOrder',102, ...                       % filter order
    'PassbandFrequency1',900, ...               % band edge
    'StopbandFrequency1',950, ...               % band edge
    'StopbandFrequency2',1050, ...              % band edge
    'PassbandFrequency2',1100, ...              % band edge
    'SampleRate',fs)                            % sampling rate
y = filter(bsFilter,x);

I used an FIR filter above, since FIR filters often give nicer results four audio, due to their linear phase response (i.e. all frequencies experience equal delay). I used order=52 because a sharp notch requires a high order.

If you also want a notch at 5 kHz, then create another bandstop filter, and pass the already-filtered signal through it.

[William Rose]

@Nathaniel,

Here the input is a linear chirp from 100 to 2100 Hz, from time 0 to 10 seconds. Therefore it will be at 1000 Hz at time t=4.5 seconds. Let's filter it with a 1 kHz notch filter.

% Create chirp signal
fs=1e4;                     % samping rate (Hz)
t=0:1/fs:10;                % time base for 10-second-long signal
x=chirp(t,100,10,2100);     % chirp signal
% create 100 Hz notch filter
bsFilter = designfilt('bandstopfir', ...       % response type
    'FilterOrder',102, ...                       % filter order
    'PassbandFrequency1',900, ...               % band edge
    'StopbandFrequency1',950, ...               % band edge
    'StopbandFrequency2',1050, ...              % band edge
    'PassbandFrequency2',1100, ...              % band edge
    'SampleRate',fs);                            % sampling rate
% Filter the chirp
y = filter(bsFilter,x);
% Plot original and filtered signals
figure
subplot(211), plot(t,x,'-b'); ylabel('x(t)'); grid on
subplot(212), plot(t,y,'-b'); ylabel('y(t)'); grid on; xlabel('Time (s)')

It attenuates the soignal around time t=4.5, which is when the frequency is 1 kHz. The ripple in the output reflects the ripples in the filter's frequency response.

[Nathaniel]

Passing the original time domain signal and dividing by fs in the t value ended up working. Thank you for the information.

[William Rose]

Great!