MATLAB Fourier Transform and Butterworth Filtering

[Lonari]

Hi,
I have a sound file, with noise, and I'm trying to remove the noise.
The way I've been trying to do this, is to do an FFT, and plot to find frequency, and do a low pass filter to remove the higher frequency. Here are the steps I've been using:

%Frequency f
f=fs/2*linspace(0,1,1024/2+1);
%Discrete Fourier Transform y
y=fft(data);
%Butterworth Filter with order 2, cuttoff Frequency 4500, sample frequency 11025
[B,A]=butter(2,5000/11025);
plot(f,2*abs(y(1:1024/2+1)))
filtdata=filtfilt(B,A,data);
z=fft(filtdata);
plot(f,2*abs(filtdata(1:1024/2+1)));
sound (filtdata)

It does not seem to work. I guess I can upload the sound file if needed.

 

[thylacine1975]

Hiya Lonari,

Without knowing the spectrum of the signal + noise of your sound file I'm guessing, but two things are apparent from your script:

1. You're using the Butterworth filter. While this might be OK, the Butterworth provides a maximally flat passband response by sacrificing rolloff rate. As such, it isn't the best weapon to use for separating features that are spectrally close together. An Elliptic or Tschebychev filter would probably be a better choice. The human ear is quite forgiving when it comes to tolerating passband ripple (for intelligibility, *not* hi-fi audio , so you can generously trade passband flatness for a steeper filter skirt.

2. The corner frequency of your low-pass filter is quite close to Nyquist, and thus the "filtered" region of spectrum is quite small. If you really do need the filter to suppress only a small spectral region, increase the number of poles/filter order (since group/delay/causality isn't a concern here) and see if that helps too.

MATLAB's fdatool is quite a useful aid to visualisation. I imported your current filter ...
>> [B,A]=butter(2, 5000/11025) % Fs = 11025 Hz, Fc actually = 2500 Hz (Not 5000 as you may have intended?)
and here's the response:


I have previously helped someone else with a problem similar to yours. There may be some extra clues/useful tips in the previous thread here: https://www.edaboard.com/threads/250442/
(although, for some unknown reason, they heavily edited it afterwards and deleted all of their text...?)

Good luck

 

[Lonari]

Hi thylacine! Thanks for the response, been waiting a long time for someone to respond haha!

Attached the wav file here:
**broken link removed**

when I do the fft:
f=fs/2*linspace(0,1,1024/2+1);
y=fft(data,1024);
plot(f,2*abs(y(1:1024/2+1)))

I get this:

I'm trying to get rid of that peak after 5kHz thats why I used the 5000/11025.
I'm not sure if I'm doing this right. In any case, if you have the time, I attached the wave file above maybe you can look at.
Thank you so much!

I understand that butterworth may not be the best filter to use, but since I'm trying to cutt off that huge spike, I dont think it's too bad in this case. Also if I increase the order, wouldnt butterworth be adequate anyway? Thanks again!

 

[thylacine1975]

Ah - now I can see the spectrum, low pass filtering makes sense
[I can't seem to download the .wav file though - no matter]

Aye, the Butterworth filter (with sufficiently high order and/or a low enough Fc) will do the job - the only problem is that it will also attenuate more of your desired signal spectrum than a sharper filter would (which may/may not be a problem, depending on what the desired signal spectrum is). Given there's only one narrowband tone though, I'd recommend a notch filter. Try:

filt = fdesign.notch('N,F0,BW', 10, 5000, 100, 11025); % 5 kHz notch, ~100Hz wide
filtData = filter(design(filt), data);

Let me know how it goes!

 

[Lonari]

Ah - now I can see the spectrum, low pass filtering makes sense
[I can't seem to download the .wav file though - no matter]

Aye, the Butterworth filter (with sufficiently high order and/or a low enough Fc) will do the job - the only problem is that it will also attenuate more of your desired signal spectrum than a sharper filter would (which may/may not be a problem, depending on what the desired signal spectrum is). Given there's only one narrowband tone though, I'd recommend a notch filter. Try:

filt = fdesign.notch('N,F0,BW', 10, 5000, 100, 11025); % 5 kHz notch, ~100Hz wide
filtData = filter(design(filt), data);

Let me know how it goes!

I will try that and let you know how it goes! Thanks a lot.
Sorry about the upload, but I managed to upload the file again - this time so you'd be able to download it if you'd like.

This is the unfiltered signal
**broken link removed**

Once I apply the first filter and remove that peak, I am supposed to get this:
**broken link removed**

It has several farm animals in it, and what I'm trying to do is to separate the animal sounds.
I will try the notch filter as you mentioned, but in any case, my code on my first post does not work with the carrier_v2.wav - atleast I dont know why I don't hear anything but static when I play it. I'm not sure if I'm using the butterworth wrong, or if I'm using fft in the wrong places or what haha.

Thanks again for your help!

 

[thylacine1975]

Oh, I just noticed - when you play the sound in MATLAB, don't forget to specify the sampling frequency argument:

[From the on-line help]

SOUND(Y,FS) sends the signal in vector Y (with sample frequency
FS) out to the speaker on platforms that support sound. Values in
Y are assumed to be in the range -1.0 <= y <= 1.0. Values outside
that range are clipped. Stereo sounds are played, on platforms
that support it, when Y is an N-by-2 matrix.

SOUND(Y) plays the sound at the default sample rate of 8192 Hz.

 

[Lonari]

So, I did this after taking your advice:

sound (data,fs)
f=fs/2*linspace(0,1,1024/2+1);
[B,A]=butter(2,(5000/11025));
filtdata=filtfilt(B,A,data);
sound(filtdata,fs)

I didn't do the fft for anything since I already know what frequencies I need to remove.

But even after using sound(filtdata,fs) I don't get the desired output.

I also used the notch filter but still just get static:

filt = fdesign.notch('N,F0,BW', 10, 5000, 100, 11025);
filtdata1 = filter(design(filt), data);
sound (filtdata1,fs)

 

[enjunear]

I'm not sure if I'm using the butterworth wrong, or if I'm using fft in the wrong places or what haha.

FFT and the filtering appear to be separate actions in your original file. One should not affect the other.

Here's a breakdown of your original code and what it should be doing (in pseudo-code).

%Frequency f
f=fs/2*linspace(0,1,1024/2+1); ---> creates an array of frequencies (so the FFT result is more readable)

%Discrete Fourier Transform y
y=fft(data); ----> take "data" and figure out the frequency content, save it as "y"

%Butterworth Filter with order 2, cuttoff Frequency 4500, sample frequency 11025
[B,A]=butter(2,5000/11025);  ---> generate filter coefficients for a Butterworth filter
plot(f,2*abs(y(1:1024/2+1)))  ---> plot FFT output vs. the frequency bins (shows your signal spectrum)
filtdata=filtfilt(B,A,data); ---> use butterworth filter coefficients and applies the filter to "data", saves it into "filtdata"

z=fft(filtdata);  ---> take "filtdata" and figure out the frequency content, save it as "z"
plot(f,2*abs(filtdata(1:1024/2+1))); ---> plot "filtdata" frequency content vs. frequency bins

sound (filtdata)  --> play filtered sound file (filtdata)

You could remark out (or delete) everything related to the FFT and plotting and all you'd be left with is 3 lines.

[B,A]=butter(2,5000/11025);  ---> generate filter coefficients for a Butterworth filter
filtdata=filtfilt(B,A,data);  --> apply filter defined by coefficients in A and B
sound (filtdata)  --> play filtered sound file (filtdata)

All of the FFT-related code is just to allow you to visualize what the signal looks like before and after you apply the filter. Creating the vector 'f' just gives you an easy to read plot, so that the fft data lines up with a useful number (like the large tone at 5 kHz).

 

[Lonari]

Thanks enjunear,
I do see it clearly now, the first code I posted was when I first dabbled in it and It wasn't too clear to me then. Even with the 3 lines of code that actually do do the filtering I don't get what I'm looking for. The static still seems to be there and that's the main problem I'm having. I did play the sound with the sample frequency, i.e sound (filtdata, fs) and I still get static, so I'm not sure how to move forward.

 

[thylacine1975]

I see what's wrong though... there aren't any farm animals hiding in the data of carrier_v2.wav!

Behold the results of:
>> figure; pwelch(data, [], [], [], fs);


...so when you filter away the single tone...
(using the correctly tuned notch filter: filt = fdesign.notch('N,F0,BW', 10, 5264, 100, 11025); % 5.264 kHz notch, ~100Hz wide)
... you are left with nothing *but* "static".

(note that it doesn't look anything like the spectrum of farm_v2.wav)

Unless there's something far more complicated like broadband masking (which isn't going to be recovered with a simple linear filter) happening here, someone's yanking your chain!
Good luck getting to the bottom of it

(P.S. The filename "carrier" also suggested to me the possibility that perhaps the farm animal sounds were carried as modulation on the 5.264 kHz carrier signal - but we can discount this possibility as a) we don't see the characteristic shape of modulation sidebands in the data, and b) there's not enough bandwidth available about the 5.264 kHz carrier in the 11025 Hz sampled data)

 

[Lonari]

Haha! Possibly you are right! I've been torturing myself all weekend trying to get this to work! I'm gonna talk to my professor in about 7hrs when I get to school in the morning, and see whats up because he might have uploaded the wrong file on there! I'll keep you updated! Thanks for all your help!

- - - Updated - - -

So I really feel silly about this, but he said that I'm supposed to use the farm_v2 and separate out the farm animal sounds! I'm not sure why he gave us the carrier_v2 to begin with!
So this is my initial fft plot for the farm_v2:

And after applying the butterworth filter I get this:

What filter would you recommend to separate out those animal sounds?
If I could maybe make use of windowing I could just specify a range and get those sounds out? Does that even work? Or should I use the butterworth at higher orders?
Now that I know my code wasn't wrong to begin with I'm a little more confident!

 

[enjunear]

I could hear three animals in that file (cow, dog, bird... and an annoying 5 kHz whistle, which should now be gone). The cow is predominantly low frequencies, the dog is somewhat in the middle of the band, and the bird is up high. I'd say you could run three filters on this and separate those three ranges. So, a low pass, bandpass and a highpass.

Look at thylacine's welch power spectral density plot.... it looks like all of the power in the signal is contained below 3 kHz. I'd say start with a 1 kHz lowpass (0-1 kHz), a 1 kHz wide bandpass centered at 1.5 kHz (1-2 kHz), and a highpass with a cutoff around 2 kHz (2 kHz+). You might swap out the highpass for another bandpass from 2-3 kHz, in case there is a lot of extra high frequency junk that comes through after you filter off the low end signals. I could be wrong and the bird sound contains some higher frequency content up past 3 kHz that you may need to keep (i.e. highpass filter).

Give that a shot and see what you can hear in each of the bands.

--------

As an afterthought, check out the time periods when each animal makes it's sound. You can also truncate the file and then run your filter on the shorter section of sound. That way to leverage both the sound separations in frequency and time.

 

[Arka Das]

Hi, Currently I am also working on the same problem i.e. filtering the farm noise and identification of animals as defined above. Hence it will be very helpful to me if anyone can kindly help me out with the coding and algorithm as I am unable to detect any of the farm animal till now.
Thanks