audio - Audioread in matlab of wav file and FFT

Question

I'm working on Matlab, I want to perform FFT on a wav file I previously recorded on Matlab as well.

fs = 44100; % Hz
t = 0:1/fs:1; % seconds
f = 600; % Hz

y = sin(2.*pi.*f.*t);

audiowrite('600freq.wav',y,fs)

This is the way I'm recording in the wav file. Now to the reading and FFT part:

[y,Fs] = audioread('600freq.wav');
sound(y)
plot(fft(y))

This is the plot of the FFT I get:

Maybe I'm missing something about the FFT, but I expected two vertical lollipops. Another thing I noticed that's wrong, is when I play the sound after reading it form the file it's longer and the pitch is significantly lower. My guess is a sampling rate problem, but I really have no idea of what to do about it.

Thanks for any help in advance.

Answer 1

That's because you're not plotting the magnitude. What you are plotting are the coefficients, but these are complex valued. Because of that, the horizontal axis is the real component and the vertical axis is the imaginary component. Also, when you use sound by itself, the default sampling frequency is 8 kHz (8192 Hz to be exact) which explains why your sound is of a lower pitch. You need to use the sampling frequency as a second argument into sound, and that's given to you by the second output of audioread.

So, try placing abs after the fft call and also use Fs into sound:

[y,Fs] = audioread('600freq.wav');
sound(y, Fs);
plot(abs(fft(y)))

Also, the above code doesn't plot the horizontal axis properly. If you want to do that, make sure you fftshift your spectra after you take the Fourier transform, then label your axis properly. If you want to determine what each horizontal value is in terms of frequency, this awesome post by Paul R does the trick: How do I obtain the frequencies of each value in an FFT?

Basically, each horizontal value in your FFT is such that:

F = i * Fs / N

i is the bin number, Fs is the sampling frequency and N is the number of points you're using for the FFT. F is the interpreted frequency of the component you're looking at.

By default, fft assumes that N is the total number of points in your array. For the one-sided FFT, i goes from 0, 1, 2, up to floor((N-1)/2) due to the Nyquist sampling theorem.

Because what you're actually doing in the code you tried to write is displaying both sides of the spectrum, that's why it's nice to centre the spectrum so that the DC frequency is located in the middle and the left side is the negative spectra and the right side is the positive spectra.

We can incorporate that into your code here:

[y,Fs] = audioread('600freq.wav');
sound(y, Fs);
F = fftshift(abs(fft(y)));
f = linspace(-Fs/2, Fs/2, numel(y)+1);
f(end) = [];    
plot(f, F);

The horizontal axis now reflects the correct frequency of each component as well as the vertical axis reflecting the magnitude of each component.

By running your audio generation code which generates a sine tone at 600 Hz, and then the above code to plot the spectra, I get this:

Note that I inserted a tool tip right at the positive side of the spectra... and it's about 600 Hz!