Class blog 2/21

ECE 3340 - Han Q. Le (c)

Help for HW 3 and classwork.

1. Linear regression of Fourier harmonics

Consider this:

Suppose we have measured various dominant frequencies from 1 to 20 and put into a paired list.
(to pair, use Transpose[{Range[Length[fList],fList}]

What we can do to find the fundamental harmonic is to do a linear regression. First, let’s plot the data:

fdataplot=ListPlot[fharmonics,PlotStyle->{Red,PointSize[0.015]}]

We want to fit a straightline through them, this is how we can do it:

fFit=LinearModelFit[fharmonics,x,x,IncludeConstantBasis->False
(* Weights-> amplitude - we don't need to use weights here *)]

ffund=fFit["BestFitParameters"][[1]]

Let’s check by plotting the fitting on top of the data:

Show[fdataplot,Plot[fFit[x],{x,fharmonics[[1,1]],Last[fharmonics][[1]]}]]

Quite a good fit, isn’t it? There is a parameter (that we will learn later on), called R-squared that indicates how good the fit is:

fFit["RSquared"]

Well, 1 is perfect, so this is quite remarkable. Don’t expect to always get this. 0.999 is quite good and that's all expected from HW.

resid=fFit["FitResiduals"]

and we can plot:

ListPlot[Transpose[{fharmonics[[   ;;,  1 ]],resid}]
, PlotStyle->{Red,PointSize[0.02]}
,Filling->Axis]

We can also plot the relative residuals:

ListPlot[Transpose[{fharmonics[[   ;;,  1 ]],resid/Abs[fharmonics[[   ;;, 2 ]]]}]
, PlotStyle->{Blue,PointSize[0.02]}
,Filling->Axis]

In other words, the deviation is <= 0.3%.

2. Find Fourier coefficients

2.1 Fourier coefficients

There are a number of algorithms with increasing sophistication and accuracy, but here, we will use the simplest and quickest: nearest maximum. 
Besides data and sampling rate, we need a list of the number of dominant harmonics, such as {2, 3, 6, 11, 8, 9, 7...}. It doesn’t have to contain every integer in a range. For example, a signal with a very weak 5th harmonic may not need the coefficient for m=5. However, if it is too bothersome having to select which one to include, which one to exclude, one can include the entire range if the added extra noises or interference are insignificant. Hence, one can just stipulate a range in low noise signals. Also, numeric order is not important. In fact, it can go from the highest to the lowest power harmonic in descending order, for example.

Go to Section 4 below to see utility  function findFourier0

Below is an example:

nharmL=Range[20];
{fouriercoef,fcpair,ftabsSpec,ftdBSpec,ftabsmx,ftdBmx}=findFourier0[{x1,srate},{nharmL,ffund}];

Proof-in-the-pudding how good the harmonic linear regression is evident in the graph below.

Show[ListPlot[ftabsSpec, PlotRange->{{0.,1.2*Max[nharmL]*ffund},All}
,PlotStyle->{PointSize[0.002],Hue[0.47,1,1]}]
,Graphics[{Hue[0.08,1,1,0.5],Thickness[0.005]
,Line[{{#[[1]],0},{#[[1]],Abs[#[[2]]]}}] & /@ fcpair}]
,Frame->True, GridLines-> Automatic
,FrameStyle-> White, GridLinesStyle-> GrayLevel[0.5]
,Epilog-> Inset[Style["dots: Fourier components
orange lines: selected coefs",18,Yellow, FontFamily->"Calibri"]
,ImageScaled[{0.75,0.85}]]
,Background-> Black, ImageSize-> 600]

We can see on dB scale:

Show[ListLinePlot[ftdBSpec, PlotRange->{{0.,1.2*Max[nharmL]*ffund},{ftdBmx-60,ftdBmx+5}}
,PlotStyle->{Thickness[0.001],Hue[0.15,1,1]}]
,Graphics[{Hue[0.85,1,1,0.5],Thickness[0.005]
,Line[{{#[[1]],ftdBmx-60},{#[[1]],20Log[10.,Abs[#[[2]]]]}}] & /@ fcpair}]
,Frame->True, GridLines-> Automatic
,FrameStyle-> White, GridLinesStyle-> GrayLevel[0.5]
,Epilog-> Inset[Style["yellow: power spectral density
vertical lines: selected coefs",18,Yellow, FontFamily->"Calibri"]
,ImageScaled[{0.75,0.85}]]
,Background-> Black, ImageSize-> 600]

2.2 Graphics you should present for HW 3- problem 1

You should have a list of frequency, and their harmonic orders in order to do linear regression. Let’s call that list fList. The graph in 2.1 above shows the spectral lines corresponding to the selected Fourier coefficients, but they are not fList. You should have a graph of fList, and this is how to do it:

{fouriercoef,fcpair,ftabsSpec,ftdBSpec,ftabsmx,ftdBmx}=
findFourier0[{x1,srate},{nharmL,ffund}];

Feel free to change the style below to your taste. As long as the graph is clear and easy to read.

ftlinplot=ListLinePlot[ftabsSpec,PlotRange-> {{0.,1.1* Max[fList]},All}
,PlotStyle->{Thickness[0.0035],Hue[0.05,1,1]}
,Filling->{1-> {0,Hue[0.08,1,1,0.4]}}
,Frame-> True, GridLines->Automatic
,FrameStyle->White, GridLinesStyle->GrayLevel[0.35]
,FrameLabel->{"frequency (Hz)","Fourier component magn."}
,LabelStyle->{16, FontFamily->"Calibri"}];
fListlingr=Graphics[{Thickness[0.001],Hue[0.45,1,1,0.8]
,Table[Line[{{f,0},{f,ftabsmx}}],{f,fList}]}];
Show[ftlinplot,fListlingr,ImageSize->500,Background-> Black]

ftdBplot=ListLinePlot[ftdBSpec
,PlotRange-> {{0.,1.1* Max[fList]},{ftdBmx-80,ftdBmx}}
,PlotStyle->{Thickness[0.0035],Hue[0.05,1,1]}
,Filling->{1-> {Bottom,Hue[0.08,1,1,0.4]}}
,Frame-> True, GridLines->Automatic
,FrameStyle->White, GridLinesStyle->GrayLevel[0.35]
,FrameLabel->{"frequency (Hz)","Fourier component (dB)"}
,LabelStyle->{16, FontFamily->"Calibri"}];
fListdBgr=Graphics[{Thickness[0.001],Hue[0.45,1,1,0.8]
,Table[Line[{{f,ftdBmx-80},{f,ftdBmx}}],{f,fList}]}];
Show[ftdBplot,fListdBgr,ImageSize-> 500,Background-> Black]

Example

fList={141.4,286.1,423.1,1842.6,1987.2,2124.2,2268.9,2405.9,2695.2,2832.2,2976.9,3121.5,3258.5,3403.2,3547.8} ;

nharmL=20; ffund=141.77;
{fouriercoef,fcpair,ftabsSpec,ftdBSpec,ftabsmx,ftdBmx}=findFourier0[{x1,srate},{mhall,ffund}];

ftlinplot=ListLinePlot[ftabsSpec,PlotRange-> {{0.,1.1* Max[fList]},All}
,PlotStyle->{Thickness[0.0035],Hue[0.05,1,1]}
,Filling->{1-> {0,Hue[0.08,1,1,0.4]}}
,Frame-> True, GridLines->Automatic
,FrameStyle->White, GridLinesStyle->GrayLevel[0.35]
,FrameLabel->{"frequency (Hz)","Fourier component magn."}
,LabelStyle->{16, FontFamily->"Calibri"}];
fListlingr=Graphics[{Thickness[0.001],Hue[0.45,1,1,0.8]
,Table[Line[{{f,0},{f,ftabsmx}}],{f,fList}]}];
Show[ftlinplot,fListlingr,ImageSize->500,Background-> Black]

ftdBplot=ListLinePlot[ftdBSpec
,PlotRange-> {{0.,1.1* Max[fList]},{ftdBmx-80,ftdBmx}}
,PlotStyle->{Thickness[0.0035],Hue[0.05,1,1]}
,Filling->{1-> {Bottom,Hue[0.08,1,1,0.4]}}
,Frame-> True, GridLines->Automatic
,FrameStyle->White, GridLinesStyle->GrayLevel[0.35]
,FrameLabel->{"frequency (Hz)","Fourier component (dB)"}
,LabelStyle->{16, FontFamily->"Calibri"}];
fListdBgr=Graphics[{Thickness[0.001],Hue[0.45,1,1,0.8]
,Table[Line[{{f,ftdBmx-80},{f,ftdBmx}}],{f,fList}]}];
Show[ftdBplot,fListdBgr,ImageSize-> 500,Background-> Black]

3. How to generate and obtain a clean spectrum

Will be discussed in class. Examples will also be shown.

4. General utility functions

4.1 Utility software for HW 3 - problem 1

findFourier0[{data_,srate_},{nharmList_,ffund_}]:=Block[{xFT,xFTAbs,xFTdB,xFTmax
,Δf,fR,xFTspec,xFTdBspec,xFTdBmax,fsel0,fsel1,ipos,coef}
,
xFT=Take[Drop[Fourier[data-Mean[data]],1],Floor[Length[data]/2]];
xFTAbs=Abs[xFT]; xFTmax=Max[xFTAbs];
xFTdB=20Log[10.,xFTAbs];xFTdBmax=Max[xFTdB];
Δf=1.*srate/Length[data];
fR=Range[Length[xFT]]*Δf;
xFTspec=Transpose[{fR,xFTAbs}];xFTdBspec=Transpose[{fR,xFTdB}];
fsel0=Table[Select[xFTspec,Abs[#[[1]]-m*ffund]≤ 2 Δf &],{m,nharmList}];
fsel1=Transpose[Flatten[Table[Take[SortBy[u,Last],-1],{u,fsel0}],1]][[1]]; 
ipos=Position[fR,#][[1,1]] & /@ fsel1;
coef=xFT[[#]]& /@ ipos;
{coef,Transpose[{nharmList*ffund,coef}] ,xFTspec,xFTdBspec
,xFTmax,xFTdBmax}
];


signalgen[{freq_,tdur_,srate_},coef_]:=Module[{nsampl,Δt,tarr
,fExpn,sigout}
,
nsampl=Round[tdur*srate];Δt=1./srate;
tarr=(Range[nsampl]-1.)*Δt;
fExpn=Table[Exp[I 2 π n  Mod[freq tarr,1]],{n,Length[coef]}];
sigout=Re[coef . fExpn];
{sigout
,ListLinePlot[Take[sigout,Round[5*srate/freq]]
,DataRange->{0,5*1./freq}
,Frame-> True, GridLines->Automatic
,FrameLabel->{"time (whichever 1/(unit of frequency))"
,"amplitude (a. u.)"}
,LabelStyle-> {18,Black}]}
];
envelopF[x_,td_,{itype_,τ_}]:=If[itype<1,
((E^(-x/τ))/(1+(E^(-(x-0.02)/0.05))))   
,(1/(1+(E^(-x/0.2))))  ((E^(-(x-td)/0.025))/(1+(E^(-(x-td)/0.025))))*(1+0.2Sin[2 π 3.5 x])] ;
notefreqdata=Import["https://pages.mtu.edu/~suits/notefreqs.html","Data"][[2]];

4.2 Generate music

Example

fouriercoef={1,0.5,0.25,0.25,0.,0.5};
notes={"C4","G4","C5"};
notedur=1/{2,2,2};

freq=Table[Select[notefreqdata
,StringContainsQ[StringReplace[#[[1]]," "-> ""],n] & ][[1,2]]
,{n,notes}];
beatdur=4; srate=44100;Δt=1/srate;
tdur=notedur*beatdur; (* actual time duration of each note *)
tau=0.5;  (* decay time for percussion/pluck instruments  *)
sigparm=Transpose[{freq,tdur,ConstantArray[srate,Length[freq]]}];
output=Flatten @ Table[su=signalgen[sigparm[[i]],fouriercoef][[1]];
tx=(Range[Length[su] ]-1)/srate;
su*envelopF[tx,sigparm[[i,2]],{1,tau}]
,{i,Length[sigparm]}];
ListPlay[output,PlayRange->All, SampleRate->srate]