3.10 Use FFT to display the power spectrum of the content of an audio wav file

Problem: download a wav file and display the frequency spectrum of the audio signal using FFT. The Matlab example was based on Matheworks tech note 1702.

3.10.1 Mathematica

Remove["Global`*"]; fname = "stereo.wav"; SetDirectory[NotebookDirectory[]]; ele = Import[fname,"Elements"]	{AudioChannels,AudioEncoding, Data,SampledSoundList, SampleRate,Sound}
(*listen to it *) sound = Import[fname,"Sound"]; EmitSound[sound] fs = Import[fname,"SampleRate"]	22050
now load the samples and do power spectrum data = Import[fname,"Data"]; {nChannel,nSamples}=Dimensions[data]	Out[115]= {2,29016}
py = Fourier[data[[1,All]], FourierParameters->{1,-1}]; nUniquePts = Ceiling[(nSamples+1)/2]	Out[117]= 14509

 

py = py[[1;;nUniquePts]]; py = Abs[py]; py = py/nSamples; py = py^2; If[ OddQ[nSamples], py[[2;;-1]]=2*py[[2;;-1]], py[[2;;-2]]=2*py[[2;;-2]] ]; f = (Range[0,nUniquePts-1] fs)/nSamples; ListPlot[Transpose[{f,py}], Joined->True, FrameLabel->{{"|H(f)|",None}, {"hz","Magnitude spectrum"}}, ImageSize->400, Frame->True, RotateLabel->False, GridLines->Automatic, GridLinesStyle->Dashed, PlotRange->{{0,12000},All}]

 

3.10.2 Matlab

clear all; close all; 
% 
%This script plots the frequency spectrum of a wave file. 
%This method of plotting the frequency spectrum is modeled after the 
%example on Mathworks website tech note 1702 
 
[data,Fs]=audioread('stereol.wav'); 
[nSamples,nChannels]=size(data); 
waveFileLength=nSamples/Fs; 
 
N = 2^(nextpow2(length(data(:,1)))); 
 
Y=fft(data(:,1),N); 
NumUniquePts = ceil((N+1)/2); 
Y = Y(1:NumUniquePts); 
P = abs(Y); 
P = P/length(data(:,1)); 
P = P.^2; 
 
if rem(N, 2) % odd nfft excludes Nyquist point 
  P(2:end) = P(2:end)*2; 
else 
  P(2:end -1) = P(2:end -1)*2; 
end 
 
f = (0:NumUniquePts-1)*Fs/N; 
plot(f,P); 
title(sprintf('Power Spectrum of a wave file')); 
xlabel('Frequency Hz'); 
ylabel('|H(f)|'); 
grid; 
print(gcf, '-dpdf', '-r600', 'images/matlab_e52_1');