Fourier Tutorial Series - 2
time-dependent signals

ECE 3340  & generic                         Han Q. Le (c)

1. Fourier synthesis for time-dependent phenomena

1.1 Time-dependent signals

We have seen this previously

coef1={0.18625,0.110645,-0.02088,0.04475
,0.01658,0.00237,-0.01098,0.003489
,0.00573,0.00473,-0.00203,0.00018,0.00027
,0.00139,-0.00005,0.00154,0.00054};

coef2={0.5118,-0.0237,-0.2368,0.1724
,-0.0215,0.0301,-0.0352,-0.0098
,0.006,0.0036,0.014,-0.0092,-0.0008
,-0.0068,0.0025,0.0055,-0.0009} ;

Both are made up of the same basis functions. Only the coefficients are different.

What if we make them to be time-signal instead of space-signal? Besides the “look,” what are other differences?

Numerical approach: coding a function using Block or Module to generate a time signal

Suppose we wish to generate a time-periodic signal, what do we need to know:
1- the frequency: call it freq
2- the amplitude and phase: this is represented as the complex amplitude of a phasor:
                        A (e^(i 2 π f t)) = |A| (e^(i ϕ))  (e^(i 2 π f t))
3- the sampling rate: how often do we sample the signal as output? t is continuous, but there is no such a thing as a continuous numerical array, we have to sample time with a finite set of discrete sampling point. Let’s call it srate
4- duration: how long in time is the signal? The signal may be forever, but numerically, there is a finite size of an array of numbers to be generated.

So, below is an example of such a function, we call it signalgen for signal generator:

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}]}
];

Example how to use it: we will generate a random signal to test.

nharmonics=RandomInteger[{4,16}];
coefLogamp=RandomReal[{-1,0},nharmonics];
coefϕ=π  RandomReal[{-1,1},nharmonics];
coef=(E^(coefLogamp+ I coefϕ));
freq=440.;tdur=1.;srate=22050;
u=signalgen[{freq,tdur,srate},coef];
u[[2]]

Now, we can apply to the two sets of coefficients:

srate=22050; freq=440. ; tdur=1. ;
sigparm={freq,tdur,srate};
coef1={0.18625,0.110645,-0.02088,0.04475
,0.01658,0.00237,-0.01098,0.003489
,0.00573,0.00473,-0.00203,0.00018
,0.00027,0.00139,-0.00005,0.00154,0.00054};
coef2={0.5118,-0.0237,-0.2368,0.1724
,-0.0215,0.0301,-0.0352,-0.0098
,0.006,0.0036,0.014,-0.0092,-0.0008
,-0.0068,0.0025,0.0055,-0.0009} ;
coef1S=Drop[coef1,1]; coef2S=Drop[coef2,1];
s1=signalgen[sigparm,coef1S];
s2=signalgen[sigparm,coef2S];
Grid[{Show[#,ImageSize->350,AspectRatio-> 0.3]&
/@{s1[[2]],s2[[2]]}}]

1.2 Test signal difference as sound waves

Grid[{Show[#,ImageSize->350,AspectRatio-> 0.3]&
/@{s1[[2]],s2[[2]]}
,ListPlay[#,PlayRange-> All,SampleRate->srate] & /@ {s1[[1]],s2[[1]]}
}]

Both have the same frequency, why do they sound difference?

1.3 Adding some more features

envelopF[x_,td_]:=(1/(1+(E^(-x/0.2))))  ((E^(-(x-td)/0.025))/(1+(E^(-(x-td)/0.025))))*(1+0.2Sin[2 π 3.5 x]) ;

freq={261.63,392.,523.25};  tdur={0.75,0.75,0.75} ;
srate=44100; 
sigparm=Transpose[{freq,tdur,ConstantArray[srate,3]}];

coef1S={0.110645,-0.02088,0.04475
,0.01658,0.00237,-0.01098,0.003489
,0.00573,0.00473,-0.00203,0.00018
,0.00027,0.00139,-0.00005,0.00154,0.00054};
coef2S={-0.0237,-0.2368,0.1724
,-0.0215,0.0301,-0.0352,-0.0098
,0.006,0.0036,0.014,-0.0092,-0.0008
,-0.0068,0.0025,0.0055,-0.0009} ;
s1L=Flatten @ Table[su=signalgen[sigparm[[i]],coef1S][[1]];
tx=(Range[Length[su] ]-1)/srate;
su*envelopF[tx,sigparm[[i,2]]]
,{i,Length[sigparm]}];

s2L=Flatten @ Table[su=signalgen[sigparm[[i]],coef2S][[1]];
tx=(Range[Length[su] ]-1)/srate;
su*envelopF[tx,sigparm[[i,2]]]
,{i,Length[sigparm]}];

Grid[{Show[#,ImageSize->350,AspectRatio-> 0.3]&
/@{s1[[2]],s2[[2]]}
,ListPlay[#,PlayRange-> All,SampleRate->srate] & /@ {s1L,s2L}
}]

2. Review - the nature of sound

1.1 What is the nature of sound?

Sound in air

Sound in solid (on a solid surface)

See Tutorial 2.5

1.2 We made sound previously in Tutorial 2 (using computer)

Given:

notes={"D5","D5","E5","D5","G5","F#5"
,"D5","D5","E5","D5","A5","G5"};
notedur=1/{8,16,4,4,4,2
,8,16,4,4,4,2};
tau=0.4;  beatdur=2. ; srate=11025;
tdur=notedur*beatdur; (* actual time duration of each note *)

Answer

notes={"D5","D5","E5","D5","G5","F#5"
,"D5","D5","E5","D5","A5","G5"};
notedur=1/{8,16,4,4,4,2
,8,16,4,4,4,2};
notefreqdata=Import["https://pages.mtu.edu/~suits/notefreqs.html","Data"][[2]];
freq=Table[Select[notefreqdata
,StringContainsQ[StringReplace[#[[1]]," "-> ""],n] & ][[1,2]]
,{n,notes}];
beatdur=2. ; srate=11025;Δt=1/srate;
tdur=notedur*beatdur; (* actual time duration of each note *)
tau=0.4; 
output=Flatten @ Table[Table[(E^(-t/tau)) Sin[2 π freq[[i]] t],{t,0.,tdur[[i]],Δt}],{i,Length[notes]}];
ListPlay[output, SampleRate->srate]

Export["HBtoyou.mp3",]

The takeaways: How to get and look up reference info, data.

1.2.1 Use function Map.

We can define a function to generate a sound signal. Then use Map to play from a list.

sndP[parm_]:=Table[(E^(-t/parm[[3]])) Sin[2 π parm[[1]]  t],{t,0.,parm[[2]],parm[[4]]}] ;
pdat=Transpose[{freq,tdur,ConstantArray[tau,Length[freq]],ConstantArray[Δt,Length[freq]]}];
outP=Flatten @ Map[sndP,pdat];  (* Map syntax is identical to Python  *)
ListPlay[outP, SampleRate->srate]

sndV[parm_]:=Table[(((1-(E^(-3*t/parm[[3]])))/(1+(E^(5*(t-0.9parm[[2]])/parm[[3]])))) ) 
(Sin[2 π parm[[1]]  t]
+0.2 Sin[2 π 0.5*parm[[1]]  t]
+0.4Sin[2 π 2 parm[[1]]  t]
+0.3Sin[2 π 3 parm[[1]]  t]),{t,0.,parm[[2]],parm[[4]]}] ;

pdat=Transpose[{freq,tdur,ConstantArray[tau,Length[freq]],ConstantArray[Δt,Length[freq]]}];outV=Flatten @ Map[sndV,pdat];
ListPlay[outV, SampleRate->srate]

Export["HBtoyou2.mp3",]

3. Experimental sound analysis - harmonics

For this work, we need to import this app:

http://courses.egr.uh.edu/ECE/ECE3340/Apps/ECEgen_APP_PSD_analysis_M12_v0.nb

Notes about sound data analysis

import sound files from previous or elsewhere recordings

Example:

audiodataGet[xaudio_]:={AudioData[#][[1]],QuantityMagnitude @ AudioSampleRate[#] }&  @  xaudio;

dirnm= (* Directory of the files to be imported *) ;
SetDirectory[dirnm];  
filenm={"filename1","filename1","etc",,, };
recimp=Import[#<>".mp3"] & /@ filenm ; (* many other formats besides mp3 are accepted *)
(* wav (for Windows) and mp3 are two most consistently reliable formats that never had reported problems
in previous classes  *) 
{snddata,srate}=Transpose[audiodataGet /@ recimp];
If[Length[DeleteDuplicates[srate]]>1,"warning: different sample rates"
,srate=srate[[1]]];
(* if you import just one file, use this  *)
{whateverdata,srateofit}=audiodataGet  @ Import["whateverfile.filetype"];

entering data into the App

See instruction in class: Click “Direct data”, enter {dataname, samplerate} then, DO NOT HIT ENTER, just click elsewhere on side of the entry box of the app.

If analyzing sound data, make sure click this    to select sound:

After toggling, the icon for sound   should appear:

Follow instruction in class for analyzing the signal.

Do not open multiple windows of the same function, for example, do not open multiple linear scale spectrum windows. Mathematica computes dynamically (like games) and it would treat two windows sharing the same variables as a game between them.

The above applies to all cases in general , not just here.

3.1 Experiment 1: whistle

3.2 Experiment 2: a constant pitch vocal sound

Example

3.3 Experiment 3: same pitch, different vowels

from above, example 1:

another sound

different sound

and so on ...

4. Experimental sound analysis - musical instruments

4.1 Demo - spectrogram

HW: analyze instrument spectral difference.

instrument= {"Violin","AltoSax","Piano","Shamisen"};
soundtest=Sound[SoundNote["A",1,# ]] & /@ instrument

soundtest={,,,};

instrumdat=AudioData /@ soundtest;
(* fname=# <>".csv" & /@ instrument ; *)
(* Export[#[[1]],#[[2,1]]] & /@ Transpose[{fname,xdat}] *)

vocal= {"VoiceAahs","VoiceOohs"};
sndchoir=Sound[SoundNote["A",1,# ]] & /@ vocal

sndchoir={,};
choirdat= AudioData[#][[1]] & /@ sndchoir;

Spectrogram[#[[1]]] & /@ instrumdat

{,
,
,}

Spectrogram[#] & /@ choirdat

4.2 Analysis of imported or pre-generated data (HW)

soundtestInst={,
,,};

xdat=AudioData /@ soundtestInst;
Dimensions /@ xdat

Use Spectral analysis app.

5. Experimental sound analysis - non-harmonic

5.1 Experiment 5.1: signal in noise

Search for needles in the haystack by doing this:

srate=44100; tarr=Range[0,1.,1/srate]; 
fs={523.25,659.25,783.99};
fexp=Table[(E^(I 2 π f tarr)),{f,fs}];
amp=RandomReal[{-0.1,0.1},Length[fs]];
sigsmall=Re[ amp . fexp];
noise=RandomVariate[NormalDistribution[0.,1],Length[sigsmall]];
signoise=sigsmall+noise ;
Grid[{{"signal","noise","signal+noise"}
,ListPlay[# ,PlayRange->All, SampleRate->srate] & /@ {sigsmall,noise,signoise}}]

Use the app to get the power spectral density of the final signal. On the linear scale, we see at least two definitive peaks above the noise floor.

Check their amplitude:

amp

Compare them with noise:

hmax=(1.2/(1.*(Sqrt[2 π])));
Show[Histogram[noise,Automatic,"PDF"]
,Graphics[{Thickness[0.01]
,Table[{Hue[((i-1)/(Length[fs]+1)),1,1,0.7],Line[{{amp[[i]],0},{amp[[i]],hmax}}]},{i,Length[amp]}]}]
]

In other words, the signal amplitude is so small: {-0.033564,-0.093344,0.0016925} that they look just like zero compared to noise. Yet, the spectrum function found them.

5.2 Experiment 5.2: random, ambience sounds

SystemDialogInput["RecordSound"]

ambsn=AudioData[];
Dimensions[ambsn]

ListPlay[ambsn[[1]]]

5.3 Experiment 5.3: the sound of silence

Get into system control of sound, set the microphone level to zero:

sn=SystemDialogInput["RecordSound"]

xsnd=AudioData[sn];
ListPlay[xsnd[[1]],PlayRange->All,SampleRate->22050]

Why do we hear the sound of silence?

6. Experimental sound synthesis/analysis - random noise

πsound=Flatten[Drop[RealDigits[N[Pi,22050]],-1]];
ListPlay[πsound,SampleRate->22050]

esound=Flatten[Drop[RealDigits[N[E,22050]],-1]];
ListPlay[esound,SampleRate->22050]

randmsnd=RandomReal[{0,1},22050];
ListPlay[randmsnd,SampleRate->22050]

Use app to analyze spectrum

Wait, why there is no frequency? If there is no frequency, why is there a plot with ALL frequencies? what is the nature of the plot? where does the plot come from?

Lecture break - concept of spectrum