フィルターのクリア

Passing audio file through a notch filter transfer function

11 ビュー (過去 30 日間)
Nathaniel
Nathaniel 2024 年 3 月 30 日
コメント済み: William Rose 2024 年 3 月 30 日
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.

採用された回答

William Rose
William Rose 2024 年 3 月 30 日
編集済み: William Rose 2024 年 3 月 30 日
[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.
  3 件のコメント
Nathaniel
Nathaniel 2024 年 3 月 30 日
Passing the original time domain signal and dividing by fs in the t value ended up working. Thank you for the information.
William Rose
William Rose 2024 年 3 月 30 日
Great!

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeDigital Filter Analysis についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by