Using Discrete Fourier Transforms

APAM 1601y
Columbia University

Introduction

This notebook illustrates two typical applications of the DFT/FFT: filtering (or convolution) and creating spectrograms.

Importing Sampled Data

Mathematica can import common file formats. AIFF (or .aif) stands for "audio interchange file format" is one of the most widely used format for storing and transmitting sampled sound. (Another format, that uses lossy data compression, is MP3 (or MPEG-1 Audio Layer-3). Mathematica can not import an MP3 file. MP3 generally reduces data size by a factor of 12 while retaining the perception of high fidelity.)

Set Directory Path and Audio File Import

In order to import a file, the path to the directory/folder containing the data must be specified.

Directory[]

SetDirectory["~/Desktop/"]

s = Import["Do_Not_Go_Gentle_64kb.aif"]

Show[s]

Examine Sound Object

s // Shallow
s // Dimensions

Short[ InputForm[s]]

Shallow[InputForm[s],5]

rate = Part[s,1,2](* Sample Rate: Samples Per Second *)

1/rate // EngineeringForm(* time per sample *)

Shallow[Part[s,1,1,1]] (* Data List *)

Length[Part[s,1,1,1]]

ListPlot

sData = Drop[Part[s,1,1,1],1]; (* Make exactly 2^16 = 64 kS *)

nData = Length[sData]

{Max[sData], Min[sData]}

ListPlot[Take[sData,{1,nData,10}],PlotRange->All]

Closer examination, within a narrow range of samples,…

Floor[nData{0.5,0.51}]

Take[sData, Floor[nData{0.5,0.51}]] // Shallow

ListPlot[Take[sData,Floor[nData {0.5,0.51}]],PlotRange->All,Joined->True,PlotLabel->"Digital Levels and Intervals"]

Fourier Transform

In Mathematica, Fourier[…] can be applied to sound data and used for digital signal processing.

(sf = Fourier[sData]; ) // Timing

sf // Shallow

When Fourier[…] is applied to a sequence of Real numbers, it returns a list of Complex numbers. We'll plot these functions using the special command.

ListPlotComplex[x_List,options___] := Module[{p1,p2},
p1 = ListPlot[Re[x],PlotStyle-> Blue,DisplayFunction -> Identity,options];
p2 = ListPlot[Im[x],PlotStyle-> Red,DisplayFunction -> Identity,options];
Show[p1,p2,DisplayFunction -> $DisplayFunction]]

ListPlotComplex[Take[sf, {1, Length[sf], 10}], PlotRange-> All];

ListPlotComplex[Take[sf,{1,500}],PlotRange->All,Joined->True];

ListPlot[Log[Abs[Take[sf,{1,Floor[((Length[sf])/2)],5}]]],PlotRange->{1.`,-6.`}]

frequencies = Table[i rate/nData, {i, 0, nData/2 }] ~Join~ Reverse[Table[-i rate/nData, {i, 1, nData/2-1}]] // N;

sff = Transpose[{frequencies, sf}];

ListPlot[Apply[{#1,Log[Abs[#2]]}&,Take[sff,{1,Length[sff],10}],{1}],PlotRange->All,PlotLabel->"Frequency Spectrum Hz"]

sDataf = InverseFourier[sf];
sDataf // Shallow

Max[Abs[sData - sDataf]]

Rayleigh Energy Theorem (Parseval's Theorem)

sData . sData

sf . Conjugate[sf] // Chop

Creating a Filter

The DFT transform signals sampled in time into signals sampled in frequency. A filter eliminates specific bands in frequency space. In this section, a low-pass filter and a high-pass filter are applied to Dylan Thomas's reading.

Low Pass Filter

5000 rate /nData

lowpass[j_,lp_:5000] := If[j <lp  || j > nData - lp,1.0,0.0]

lpFilter = Table[lowpass[i],{i, 1, nData}];

ListPlot[lpFilter]

sfLP = sf lpFilter;

(sfLP . Conjugate[sfLP] // Chop) / (sData .sData) (* energy fraction retained *)

ListPlotComplex[Take[sfLP, {1, Length[sf], 10}], PlotRange-> All];

Count[sfLP,c_/;Abs[c]>0]

xfLP = InverseFourier[sfLP];

xfLP // Shallow

Max[Im[xfLP]]

Max[Re[xfLP]]

ListPlay[Re[xfLP], SampleRate -> rate]

Band Pass Filter

1000 rate /nData

bandpass[j_,lp_:5000,hp_:1000] := If[(hp < j < lp) || (nData - lp < j < nData - hp), 1.0, 0.0]

bpFilter = Table[bandpass[i],{i, 1, nData}];

ListPlot[bpFilter]

sfBP = sf bpFilter;

(sfBP . Conjugate[sfBP] // Chop) / (sData .sData) (* energy fraction retained *)

ListPlotComplex[Take[sfBP, {1, Length[sf], 10}], PlotRange-> All];

Count[sfBP,c_/;Abs[c]>0]

xfBP = InverseFourier[sfBP];

ListPlay[Re[xfBP], SampleRate -> rate]

Time-Domain Filters

A filter can also be applied using ListConvolve[…]. In general, convolution requires much more computation than the FFT. Indeed, many convolution routines use the FFT to speed execution. In this section, the "time-equivalent" signals of the filters used above are computed using the IDFT.

fLPtime = InverseFourier[lpFilter];

ListPlot[Re[fLPtime],PlotRange->All,PlotLabel->"Low Pass Filter (time)"]

ListPlot[Take[Re[fLPtime],{1,100}],PlotRange->All,PlotLabel->"Low Pass Filter (time)"]

Plus @@ fLPtime

fBPtime = InverseFourier[bpFilter];

ListPlot[Re[fBPtime],PlotRange->All,PlotLabel->"Band Pass Filter (time)"]

ListPlot[Take[Re[fBPtime],{1,100}],PlotRange->All,PlotLabel->"Band Pass Filter (time)"]

Plus @@ fBPtime

ListConvolve

?ListConvolve

lista = Table[a[i], {i, 1, 10}]

ListConvolve[{1, -2, 1}, lista, 2]

lpKernel = Table[If[i<=8, i-1, 16 - i], {i, 1, 16}]/56

Plus @@ lpKernel

ListPlot[lpKernel]

skLP = ListConvolve[lpKernel,sData,8];

ListPlay[Re[skLP], SampleRate -> rate]

pk1 = ListPlot[Log[Abs[Take[sf, {1, Floor[Length[sf]/2], 5}]]], PlotRange-> {1.0, -9.0}, PlotStyle-> RGBColor[0,1,0], DisplayFunction -> Identity];

pk2 = ListPlot[Log[Abs[Take[Fourier[skLP], {1, Floor[Length[sf]/2], 5}]]], PlotRange-> {1.0, -9.0}, DisplayFunction -> Identity];

Show[pk1,pk2,DisplayFunction->$DisplayFunction,PlotLabel->"Low-Pass Convolution"]

Frequency Spectrograms

A frequency spectrogram is made by applying the DFT to a small time "window" that moves along the sequency of samples. A frequency spectrogram is also called a "time-Fourier-transform", or TFD.

tfd[x_List,windowSize_:200, padSize_:1] := Module[{i,j,n,pad, window, xw, xf, result},
n = Length[x];
pad = Table[0, {i, 1, windowSize padSize}];
window = Table[If[i ≤ windowSize/2, (i-1) ,(windowSize - i)], {i, 1, windowSize}]/(0.5 windowSize);
result = {};
Do[xw = Take[x, {j,j+windowSize-1}]; 
xw = xw - (Plus @@ xw)/n; (* remove average *)
xw = Join[pad , (xw window), pad];
xf = Fourier[xw];
result = Append[result, Abs[Drop[xf,-Floor[Length[xf]/2]]]],
{j, 1, n-windowSize -1, windowSize/2}];
Transpose[result]]

sTFD = tfd[sData];

ListDensityPlot[sTFD,Mesh->False, PerformanceGoal -> "Speed", PlotRange -> All]

ListPlay[sData, SampleRate -> rate]

windowSize=200;
padSize=1;
dataSize=windowSize (1+2 padSize);frequenciesTFD=N[Table[((i rate)/dataSize),{i,0,(dataSize/2)}]];ListPlot[frequenciesTFD,PlotLabel->"Frequency (Hz) versus Spectral Index"]

Summary

The digital Fourier transform (DFT) has become a work-horse of modern digital signal processing. Using a short sequence of sampled voice-data, we've illustrated how to manipulate and analyze the spectral content of time-varying signals.