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.

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

[Roberto Santana]

I forgot to add this function that is part of my code to make the brickwall low pass filter.

h = 2*3000*SA(2*pi*3000*t);

[Roberto Santana]

The audio file is approximately 4 seconds long and consists of an uncompressed mono-audio wav file sampled at 44.1ksps.

Answer 1

[William Rose]

@Roberto Santana,

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

[Roberto Santana]

Thank you very much for your help

[William Rose]

@Roberto Santana, you're welcome.

Answer 2

[William Rose]

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

[Roberto Santana]

I attached the audio, what I want to do with the sound file is to cut off the 8kHz tone and just have left the person talking.

[William Rose]

@Roberto Santana,

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