how filter the noise frequency of an audio file with notch filter

4 次查看（过去 30 天）
Daniel Niu
Daniel Niu 2022-11-10
评论： Star Strider 2023-5-8
Dear friend,
I want to filter the 464 Hz noise in an audio file.
but using the filtfilt functio I get the NaN.
what is my problem? Unfortunely, the community website don't support wave file.
clear;
[y,Fs] = audioread('q.wav');
signal=y;
sound(y,Fs);
pause(1)
L = length(y); % Signal length
t = 0:1/Fs:(L-1)*1/Fs; % Time vector
t=t';
X=y;
n = 2^nextpow2(L);
Y = fft(X,n);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
semilogx(f,20*log10(squeeze(P1)),'b')
hold on;
title('Magnitude spectrum of sound')
xlabel('Frequecny [Hz]')
ylabel('PSD [dB]')
R = 1e3; % resistor value [Ohms]
C = 2*171.368563054e-9; % Capacitor value [Farads]
numerator = [(R*C)^2 0 1];
denominator = [(R*C)^2,4*R*C,1];
H = tf(numerator,denominator);
[mag,phase,wout] = bode(H);
semilogx(wout(:,1)/(2*pi), 20*log10(squeeze(mag))/5, '-g','LineWidth',2); zoom on; grid on;
title('magnitude'); xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
grid on;
f0=1/(4*pi*R*C); %Hz
signal_filtered = filtfilt(numerator,denominator,signal);
hold on;
semilogx(f,20*log10(squeeze(signal_filtered(1:L/2+1))),'r')
sound(signal_filtered, Fs);
  2 个评论
Sumit
Sumit 2023-5-8
i have to design 50hz notch filter for removing noise in ECG signal which is into .wav file. how can i design it? can i use same code with change in frequency also suggest me what chnages should have.
Star Strider
Star Strider 2023-5-8
@Sumit — There is an example Remove the 60 Hz Hum from a Signal in the documentation that you can tweak. There are several other options, including bandstop (use the 'impulseResponse','iir' name-value pair for best results). Designing your own filter (preferably using ellip and its friends, including the zp2sos function) is also straightforward.
Of course, it depends on what your instructions are with respect to desigining the filter, since this appears to be a homework assignment.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2022-11-10
编辑：Star Strider 2022-11-11
Ran in:
Use the zip function to create a .zip file for the .wav file and the upload that.
Beyond that, your numerator and denominator vectors are currently continuous time (s space) vectors, not sampled z space vectors. There are several ways to convert them, and I prefer the bilinear transformation.
Assuming a sampling frequency of 44100 Hz —
format shortE
Fs = 44100;
Ts = 1/Fs;
R = 1e3; % resistor value [Ohms]
C = 2*171.368563054e-9; % Capacitor value [Farads]
numerator = [(R*C)^2 0 1];
denominator = [(R*C)^2,4*R*C,1];
[numd,dend] = bilinear(numerator, denominator, Fs)
numd = 1×3
8.8325e-01 -1.7626e+00 8.8325e-01
dend = 1×3
1.0000e+00 -1.7626e+00 7.6651e-01
f = logspace(-3, +5, fix(Fs/10));
w = 2*pi*f;
figure
freqs(numerator, denominator, w)
sgtitle('Continuous Time')
figure
freqz(numd, dend, 2^16, Fs)
set(subplot(2,1,1), 'XScale','log')
set(subplot(2,1,2), 'XScale','log')
% sgtitle('Discrete Time')
EDIT — (11 Nov 2022 at 18:18)
I just now saw the ‘q.zip’ file.
Uz = unzip('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1188068/q.zip')
Uz = 1×1 cell array
{'question1_jammed.wav'}
[y,Fs] = audioread(Uz{1})
y = 84672×1
0 4.0497e-02 8.0933e-02 1.2122e-01 1.6129e-01 2.0108e-01 2.4048e-01 2.7945e-01 3.1790e-01 3.5577e-01
Fs =
44100
L = numel(y);
t = linspace(0, L-1, L).'/Fs;
figure
plot(t, y)
grid
xlabel('Time (s)')
ylabel('Amplitude')
y_filt = filtfilt(numd, dend, y);
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTy = fft([y y_filt].*hamming(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2-1)*Fn;
Iv = 1:numel(Fv);
TrFn = FTy(:,2) ./ FTy(:,1);
figure
subplot(3,1,1)
plot(Fv, abs(FTy(Iv,1))*2)
grid
xlim([0 1E+3])
title('Original')
subplot(3,1,2)
plot(Fv, abs(FTy(Iv,2))*2)
grid
xlim([0 1E+3])
title('Filtered')
subplot(3,1,3)
plot(Fv, abs(TrFn(Iv)))
grid
xlim([0 1E+3])
title('Transfer Function')
The filter definitely notches out the signal!
.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Analog Filters 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by