Plotting FFT for audio WAV file?
303 views (last 30 days)
Show older comments
Avinash Kandalam
on 4 Nov 2020
Commented: Mathieu NOE
on 9 May 2022
Dear all,
I tried to explain as clear as possible. I want to plot "Raw FFT" file for a "WAV" file. This WAV (audio) file is acquired from a microphone for a period of 1 minute. The goal is to plot frequency distribution (0 Hz - 20 kHz).
- I want to acquire raw FFT (to see if there are any signal peaks at particular frequency) throughout 1 minute. The WAV (audio) file (only 1) is atttached to this question.
- Please help me with the code and the output graph.
I tried to execute the following code (from previous answers here) and I think it is not the right way. I think what the code shows is basically amplitude vs frequency; but not a typical FFT spectrum.
Million Thanks,
Avinash.
CODE: I tried and most likely wrong. I think as said, it is just amp vs freq, which does not give me clear picture of frequencies which lies in different ranges.
[y1,fs]=audioread('myWAVaudiofile.wav');
t=linspace(0,length(y1)/fs,length(y1));
Nfft=16777216; %power of 2 and I put a huge number so there are many data points
f=linspace(0,fs,Nfft);
X1=abs(fft(y1,Nfft));
plot(f(1:Nfft/2),X1(1:Nfft/2))
xlabel('Frequency');
ylabel ('Power???');
title ('FFT Spectrum');
OUTPUT: I only zoomed into 0-30 Hz using above code and WAV file attached (ofcourse the wole spectrum is until 20000 Hz)
6 Comments
Walter Roberson
on 4 Nov 2020
fft() is inherently a function for processing periodic signals. The assumption is that the signal continues on forever. Remember that you are decomposing into the sum of phased sine or cosine signals, and those have no endpoint.
fft() by itself is therefore not useful for processing speech, as speech consists mostly of changing information.
However, fft() can still be useful -- provided that you apply it to a window of data that is wide enough to be able to examine frequencies of interest, but which is short enough to not "drag down" the changing nature of the signal.
A typical way to handle this is to use Short-Time FFT (SFFT), or Spectrograms. Spectrograms use SFFT and some graphical representation of power, to give an idea of how power at different frequencies changes over time.
Accepted Answer
Mathieu NOE
on 4 Nov 2020
dear friends, here my little contribution to wav file spectral analysis...
enjoy !
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 8192; %
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);
% spectrogram dB scale
spectrogram_dB_scale = 100; % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums
% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [signal,Fs]=wavread('myWAVaudiofile.wav'); %(older matlab)
% or
[data,Fs]=audioread('myWAVaudiofile.wav'); %(newer matlab)
channel = 1;
signal = data(:,channel);
samples = length(signal);
dt = 1/Fs;
t = (0:dt:(samples-1)*dt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sensor_spectrum, freq] = pwelch(signal,w,NOVERLAP,NFFT,Fs);
% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);
% apply A weigthing if needed
if option_w == 1
pondA_dB = pondA_function(freq);
sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
my_ylabel = ('Amplitude (dB (A))');
else
my_ylabel = ('Amplitude (dB (L))');
end
figure(1),semilogx(freq,sensor_spectrum_dB);grid
title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sg,fsg,tsg] = specgram(signal,NFFT,Fs,w,NOVERLAP);
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg)); % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% apply A weigthing if needed
if option_w == 1
pondA_dB = pondA_function(fsg);
sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
my_title = ('Spectrogram (dB (A))');
else
my_title = ('Spectrogram (dB (L))');
end
% saturation of the dB range :
% saturation_dB = 60; % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
function pondA_dB = pondA_function(f)
% dB (A) weighting curve
n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
pondA = n./r;
pondA_dB = 20*log10(pondA(:));
end
16 Comments
Benjamin Colbert
on 7 May 2022
Mathieu,
I am trying to run your code and am running into an error. See:
https://www.mathworks.com/matlabcentral/answers/1713980-running-into-error-when-using-pwelch-for-fft
Any help would be appreciated!
Mathieu NOE
on 9 May 2022
hello Benjamin
the issue has been found - see Jonas answer
More Answers (0)
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!