How can I create a 3kHz brickwall filter to filter "movieaudio.wav" (human voice audio) that has a corrupting 8kHz tone. Comments on the code explain what I am trying to do.

24 views (last 30 days)
clear all;
%Read in the wav file and the sample rate
[x, Fs] = audioread('movieaudio.wav');
%Display the sample rate on the command line. We need this for the
%convolution to work. The integrand will use 1/Fs
Fs;
%Play the sound out the computer speaker
sound(x,Fs)
function y = SA(x)
% Trying to Create a Brickwall Filter but can't figure it out
% File: SA.M
% CALL: y = SA(x)
% This function computes sin(x)/x.
% x is assumed to be a "scalar" or a vector;
x = x(:);
y = zeros(length(x),1);
for i = 1:length(x)
if (x(i) == 0)
y(i) = 1;
else
y(i) = sin(x(i))/x(i);
end
end
% After filter is created I need to do convolution
  2 Comments
Roberto Santana
Roberto Santana on 7 Apr 2022
The audio file is approximately 4 seconds long and consists of an uncompressed mono-audio wav file sampled at 44.1ksps.

Sign in to comment.

Accepted Answer

William Rose
William Rose on 8 Apr 2022
Edited: William Rose on 8 Apr 2022
[edit: correct spelling errors]
I understand that you were trying to create an approximation of an ideal lowpass filter by using the sin(x)/x function, also known as as sinc function. Then you wanted to convolve the sinc function with the signal. That is not the easiest way to do a filter in Matlab, but it will work, and it nicely illustrates how the process works. One key is to get the scaling of the sinc function correct, both horizontally and vertically. See details below.
fs=44100; %sampling rate (Hz)
T=4; %signal duration (s)
fco=8000; %desired cutoff frequency (Hz)
Tw=0.01; %filter half-width (s)
N=fs*T; %signal duration (samples)
Nw=fs*Tw; %filter half-width (samples)
t=(0:N-1)/fs; %time base for signals (s)
tfilt=(-Nw:1:+Nw)/fs; %time base for filter (s)
h=sin(2*pi*fco*tfilt)./tfilt; %sinc function with desired cutoff freq
h(Nw+1)=2*pi*fco; %fix the NaN value in the center
h=h/(sum(h)); %normalize the sinc function, so it has unity gain at dc
%Plot the filter's impulse response.
figure;
plot(tfilt,h,'-r'); grid on;
titlestr=sprintf('Impulse response of low pass filter, f_c_o=%.0f Hz',fco);
title(titlestr); xlabel('Time (s)'); ylabel('h');
x=rand(1,N)-0.5; %create random noise signal
y=conv(x,h,'same'); %filter the signal by convolution with the filter
%Compute power spectra of raw & filtered signals.
%Uses pwelch() from Signal Proc Toolbox.
window=2048; %window width for spectrum estimate
[Pxx,f]=pwelch(x,window,window/2,window,fs);
[Pyy,f]=pwelch(y,window,window/2,window,fs);
%Plot results.
figure;
subplot(211), plot(t,x,'-r',t,y,'-b');
xlabel('Time (s)'); ylabel('Amplitude'); grid on
title('Time Domain: x(t), y(t)'); legend('x(t)','y(t)');
xlim([1 1.01]); %plot a 10 msec segment only, for clarity
subplot(212), plot(f,Pxx,'-r',f,Pyy,'-b');
xlabel('Frequency (Hz)'); ylabel('Power'); grid on
legend('|X(f)|^2','|Y(f)|^2');
title('Frequency Domain: Power spectra of x, y')
Try the code above. Good luck with your work.

More Answers (1)

William Rose
William Rose on 8 Apr 2022
@Roberto Santana, please post a sample .wav file with the 8 kHz tone. The file size will not be too large since it is just 4-5 seconds at 44.1 kKz. Thank you.
You say you want a brick wall filter but you might not like the actual results. Filters that are sharp in the frequency domain have undesirable time-domain properties, and vice versa. The uncertainty principle and all that.
You said the file is human voice. What do you want dto do with the sound file? If your goal is for someone listening to the file to understand the speech, then you may filter out frequencies below 300 Hz and above 3.4 kHz, without affecting the intelligiibility of the speech significantly.* That's how the phone system has worked for decades, including when it was all analog. A lowpass filter with a 3.4 kHz cutoff would attenuate the 8 kHz tone, and it would not take a brick wall.
*For listeners with normal hearing. Listeners with cohlear implants understand better when the speech signal includes higher frequencies.
  2 Comments
William Rose
William Rose on 8 Apr 2022
[x,Fs]=audioread('movieaudio.wav'); %read file
[b,a]=butter(4,3400/(Fs/2)); %design filter: Butterworth lowpass
%4th order, 3400 Hz cutoff
y=filter(b,a,x); %filter signal
audiowrite('movieaudio1.wav',y,Fs); %write file
This sounds a lot better.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by