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

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);

[Sumit]

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]

@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.

Answer 1

[Star Strider]

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!

.