Plotting FFT for audio WAV file?

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

1. 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.
2. 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)

[Rik]

Your question is not urgent, the reason for that is explained here.

[Avinash Kandalam]

@Walter Roberson @Rik, I see your point, thank you for sharing, sorry, my bad to push I re-updated :)

[Walter Roberson]

I suspect that you do not want fft by itself and instead want spectrogram()

[Avinash Kandalam]

@Walter Roberson: Thanks for response. 1. The code which I used was freq vs amp? 2. Spectogram, meaning freq vs dBV? I'm confused as there are many possibilites how to represent a "WAV" file in FFT.. Is there a general/basic way which people represent WAV file?

I would like to know the difference-seems fundamental understnading issue between types of FFT's with respect to MATLAB coding, could you enlighten little more in detail about what is basic code/results people do under such cicumstances. Thanks.

Another code which I found is: I used the same WAV file as attached to the original question

 [y,Fs] = audioread('myWAVaudiofile.wav');            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

X = S + 2*randn(size(t));

plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')


Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

plot(fs,P1) 
title('FFT spectrum')
xlabel('f (Hz)')
ylabel('|P1(f)|')

OUTPUT: Obviously, different output

[Walter Roberson]

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.

[Avinash Kandalam]

@Walter: Thank you for the reply. The content is basically over the signal recorded, there are different sounds happening at different frequencies.

For example, at 30 Hz and 40 Hz, if I see a strong peak (in FFT) then I believe the frequencies produced are mostly at 30 & 50 Hz. Something like this, I want to see my WAV signal to be distinguished. Thanks

Answer 1

[Mathieu NOE]

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

[Avinash Kandalam]

@Noe: Thank you for the code. I'm getting an error at hanning (line 6, in your code), I tried to use w = hann(NFFT), but still an error. Did you already run my WAV file with your code? Could you please let me know about this error? Nothing about specific error reason mentioned at "command window". Thanks

w = hanning(NFFT); 

[Avinash Kandalam]

@Noe: Just in case, I'm using "MATLAB R2020a" licence version. Not sure, if this is the issue

[Mathieu NOE]

hello

the hanning / hann window function is part of the Signal Processing Tbx

do you have it ? (check by typing ver in the command window)

if not , this is a bit of trouble because you will also not have access to pwelch and specgram functions

otherwise, the hann function is quite simple - see below :

function w = my_hann_window(N)

if ~mod(N,2)
    % Even length window
    half = N/2;
    w = calc_window(half,N);
    w = [w; w(end:-1:1)];
else
    % Odd length window
    half = (N+1)/2;
    w = calc_window(half,N);
    w = [w; w(end-1:-1:1)];
end

%---------------------------------------------------------------------
    function w = calc_window(m,n)
    %CALC_WINDOW   Calculate the generalized cosine window samples.

    x = (0:m-1)'/(n-1);
            % Hann window
            w = 0.5 - 0.5*cos(2*pi*x);

    end
end

[Avinash Kandalam]

@Noe: Thank you. Yes, the tool was missing and I installed and run it. It worked (results shown below, I just increased the Nfft value and zoomed in between 0-1000 Hz).

Question: Is there a way to show "frequency (x-axis)" and "Power- millivolts (y-axis)"- for the same WAV file? Or these are all the basic information we can obtain from "WAV file"?

I saw in "signal processing tool" there are many options to filter the signal for example, which is more tuning the WAV file. Seeing the results, it seems not very clear to understand the signal, therefore I would like to know if there are any freq vs *** plots, possible to do. Thanks a lot :)

[Mathieu NOE]

hello

my first comment : a wav file has no enginnering unit. It scaled between +/- 1 ; it's up to you to conver that in whatever engineering units .. mV or anything else.

now, between a amplitude spectrum (what I did) and a power spectrum , there is little difference , especially hen the dB scale is used

for amplitude spectrum we usually show 20*log10(fft_amplitude)

to get from amplitude to power , you have to take the squared amplitude of the fft

and in dB it is power spectrum = 10*log10(fft_amplitude ^2) = 20*log10(fft_amplitude) !

If you need more technical stuff on that, you can have a look on this (and also read the attachement)

https://fr.mathworks.com/matlabcentral/fileexchange/40002-psd-power-spectral-density-and-amplitude-spectrum-with-adjusted-fft

[Avinash Kandalam]

@Noe: Thanks a lot for your time and sharing the docs. I will work on it. Great that you supported, got little clarified. Stay safe. Bye

[Mathieu NOE]

my pleasure !

[ZHU z]

Hello, do you know how to continue to convert this audio to frequency (HZ) on the X-axis and sound pressure level (dBA/dBC) on the Y-axis? Thanks!

[Mathieu NOE]

hello @ZHU z

1/ just a general comment first - it's better to start a new post rather than mixing with the old one , this brings confusion

2/ back to your question, I am not sure if you ask how to generate the graphs you posted or what ? do you have already a code + data ?

[ZHU z]

Sorry, forgive me that I didn't make my question clear .Thank you for your reply~

My question is:

How can I draw a graph as shown in the figure based on a piece of audio (.wav)? The vertical coordinate is the sound pressure level (SPL) in A-weighting or C-weighting. I am in the same situation as @Avinash Kandalam, with only this code and a piece of audio (wav).

[y1,fs]=audioread('SAODIJI.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');

I've been troubled by this problem for a long time, can you help me? Thank you very much~

[Rik]

Please post your question in a separate thread. Feel free to post a link here to your question.

[Mathieu NOE]

hello again

@ZHU z

try this

the code has evolved (vs the one first answered) with more option but you still have the avg fft and spectrogram displays. Now you can also apply filters if needed

clc
clearvars

option_notch = 0; % 0 = without notch filter , 1 = with notch filter
    fc_notch = 50;       % notch freq

option_LPF = 0; % 0 = without low pass filter  , 1 = with low pass filter
    fc_lpf = 1000;       % LPF cut off freq
    N_lpf = 5;           % filter order

option_HPF = 0; % 0 = without high pass filter  , 1 = with high pass filter
    fc_hpf = 1000;       % HPF cut off freq
    N_hpf = 5;           % filter order

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[signal,Fs] = audioread('SAODIJI.wav');

[samples,channels] = size(signal);

% create time vector
dt = 1/Fs;
time = (0:samples-1)*dt;

%% notch filter section %%%%%%
% y(n)=G*[x(n)-2*cos(w0)*x(n-1)+x(n-2)]+[2*p cos(w0)*y(n-1)-p^2 y(n-2)]
% this difference equation can be converted to IIR filter numerator /
% denominator
if option_notch ~= 0
    w0 = 2*pi*fc_notch/Fs;
    p = 0.995;

    % digital notch (IIR)
    num1z=[1 -2*cos(w0) 1];
    den1z=[1 -2*p*cos(w0) p^2];

    % now let's filter the signal 
    signal = filter(num1z,den1z,signal);
end

%% low pass filter section %%%%%%
if option_LPF ~= 0
    w0_lpf = 2*fc_lpf/Fs;
    % digital notch (IIR)
    [b,a] = butter(N_lpf,w0_lpf);
    % now let's filter the signal 
    signal = filter(b,a,signal);
end

%% high pass filter section %%%%%%
if option_HPF ~= 0
    w0_hpf = 2*fc_hpf/Fs;
    % digital notch (IIR)
    [b,a] = butter(N_hpf,w0_hpf,'high');
    % now let's filter the signal 
    signal = filter(b,a,signal);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 2048;    % 
OVERLAP = 0.75;

% spectrogram dB scale
spectrogram_dB_scale = 80;  % 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 = 1;

%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),decim);
    Fs = Fs/decim;
    end
   signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;

%%%%%% legend structure %%%%%%%%
for ck = 1:channels
    leg_str{ck} = ['Channel ' num2str(ck) ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),plot(time,signal);grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend(leg_str);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);

% 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(2),semilogx(freq,sensor_spectrum_dB);grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
legend(leg_str);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
    [sg,fsg,tsg] = specgram(signal(:,ck),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  

    % 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+ck);
    imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
    axis('xy');colorbar('vert');grid on
    df = fsg(2)-fsg(1); % freq resolution 
    title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Frequency (Hz)');

end

%play sound
sound(signal,Fs);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

[ZHU z]

wow~Thank you so much!

I've just run the code successfully,

thanks for the great help you've given me ~

best wishes for you~

[Mathieu NOE]

tx

My pleasure !

[Benjamin Colbert]

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]

hello Benjamin

the issue has been found - see Jonas answer

Running into error when using pwelch for FFT - (mathworks.com)