FT
2022-03-18
The audio analogy with nmr
Sound waves share useful similarities to electromagnetic radiation. The are both represented as travelling sine waves. Sound waves are a periodic disturbance in physical medium (air) whereas em radiation is a periodic change in an electric field.
After a pulse excitation, the nmr signal decay contains that reflect multiple relaxations at different frequencies.
To explore the nature of wave motion and the Fourier Transform we can use audio files and R packages that create, read, and manipuate audio files (.wav).
packages:
fftw
TuneR
Below we used R to generate representations of sound. The graphs are of sound waves of different frequency, as might be detected by a microphone: the plots are amplitude of sound versus time.
library(tuneR,fftw)
Wobj1 <- 0.75*sine(880) # here is unit = freq not degree in radians
tdir <- tempdir()
tfile <- file.path(tdir, "myWave.wav")
writeWave(Wobj1, filename = tfile)
sdat1 <- tuneR::readWave(tfile)
Wobj2 <- 1.0*sine(4400) # here is unit = freq not degree in radians
tdir <- tempdir()
tfile2 <- file.path(tdir, "myWave2.wav")
writeWave(Wobj2, filename = tfile2)
sdat2 <- tuneR::readWave(tfile2)
Wobj3 <- 0.5*sine(2000) # here is unit = freq not degree in radians
tdir <- tempdir()
tfile3 <- file.path(tdir, "myWave3.wav")
writeWave(Wobj3, filename = tfile3)
sdat3 <- tuneR::readWave(tfile3)
sdat3##
## Wave Object
## Number of Samples: 44100
## Duration (seconds): 1
## Samplingrate (Hertz): 44100
## Channels (Mono/Stereo): Mono
## PCM (integer format): FALSE
## Bit (8/16/24/32/64): 32 N <- 44100
x <- seq(1,44100,1)
w1 <- sdat1@left
w2 <- sdat2@left
w3 <- sdat3@left
wt <- w1 + w2 + w3
times <- x/44100 # same as freq because t = 1 sec
length(times)## [1] 44100 # Now we set up a plotting environment of two rows and three columns (in order to hold six graphs), using par(mfrow())
par(mfrow=c(2,2))
plot(times,w1,type="l",xlab="Time/sec",ylab="Amplitude",xlim=c(0,0.01))
plot(times,w2,type="l",xlab="Time/sec",ylab="Amplitude",xlim=c(0,0.01))
plot(times,w3,type="l",xlab="Time/sec",ylab="Amplitude",xlim=c(0,0.01))
plot(times,wt,type="l",xlab="Time/sec",ylab="Amplitude",xlim=c(0,0.01))
Questions
Consider the first graph. Imagine you are the microphone, and you measure amplitude by raising and lowering your hand. Demonstrate.
According to the graph estimate how long does it take to go up and down once?
How many times will it go up and down in one second? This is the frequency.
Sound amplitude is related to compression of air. Sketch what compression of air might look like as a function of distance.
What is oscillationg up and down in the case of EM radiation?
Now we will consider how the Fourier Transform can translate a time domain representation to a frequency representation. For study of sound, and in spectroscopy, frequency representations are more useful. We use the R fftw package command FFT to do this. Now we have a graph of intensity versus frequency. We call this a transformation from the time domain to the frequency domain. For the individual time domain graphs you could easily sketch a frequency domain representation. But what about the fourth graph? Somehow the FFT is able to unscramble the waves and do the conversion.
ft1 <- abs(fftw::FFT(w1))
ft2 <- abs(fftw::FFT(w2))
ft3 <- abs(fftw::FFT(w3))
ft4 <- abs(fftw::FFT(wt))
par(mfrow=c(2,2))
plot(x,ft1/N,type="l",xlim=c(0,5000))
plot(x,ft2/N,type="l",xlim=c(0,5000))
plot(x,ft3/N,type="l",xlim=c(0,5000))
plot(x,ft4/N,type="l",xlim=c(0,5000))
The Fourier Transform
As we have seen, the FT transforms time domain data to frequency domain. This turns out to be widely important, and is key part of modern nmr (Earning Ernst the Noble Prize), other spectroscopies (FTIR), signal processing, image processing, acoustics, and even biology (“bioacoustics”).
Even though we may not encounter FTs in our math courses, we can try to get a conceptual understanding with the math we know. The key is in the idea of orthogonality.
Questions: Orthogonality
To try to understand what the FFT is doing, we can start with someting familiar, a vector represented on x y coordinates.
Draw a vector v(.75,.50) on an x y grid. Now we know that if we drop a vertical line from the end of the vector, we will find the x component. How d we find the y component? Illustrate both.
In each case, more formally, what you just did was to multiply the vector v(0.75, 0.50) by x and y unit vectors to find the scalar product.
vec1 = v(1,0) vec2 = v(0,1) vec3 <- v(2,3)
Vector multiplication
v(x1,y1) * v(x2,y2) = v(x1x2,y1y2)
- What is the scalar product of the x and y unit vectors? The unit x and y vectors are othogonal. We used orthogonality to break apart a vector into its x and y components.
Ok, what does this have to do sound waves and the fourier transform. We see that a complex waveform is composed of many individual waves. Now it turns out that the individual waves are orthogonal to each other, and we can use that fact to “pick out” the individual waves using the fourier transform.
Since the math is more complicated, let’s use R to explore. Alternatively, we can use a raphing program luki desmos.
If we multiply two functions, and find the area under the curve of the resulting function, if it is zero then those are orthogonal.
We can sometimes see if the area under the curve is zero because the positive and negative area are symmetrical.
using R or graphing calculator (desmos)
graph the product & estimate the error under the curve
sin(x)sin(x) sine(x)sin(2x) sin(x)sin(5*x)
can use estimate the area under the curve of this graph? finite or zero
The fourier transform is simply doing this over and over again, with many frequencies. when it “hits” a part of the curve that matches the frequncy (that is non-orthogonal) it has a non zero result and picks that frequency out.
[sin(x) + sin(2x) + sin(5x)] * sin(x)
or
[sin(x) + sin(2x) + sin(5x)] * sin(4x)
How do you interpret the result?
Here is example of actual sound wave & fft using R.
# **********************************************
# reading a music sample wav file
mwave = "CantinaBand3.wav" # loading the wav file
data <- tuneR::readWave(mwave)
# sampw = "sample.wav"
# sdata <- tuneR::readWave(sampw)
summary(data) # sampling rate 44,100##
## Wave Object
## Number of Samples: 66150
## Duration (seconds): 3
## Samplingrate (Hertz): 22050
## Channels (Mono/Stereo): Mono
## PCM (integer format): TRUE
## Bit (8/16/24/32/64): 16
##
## Summary statistics for channel(s):
##
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -8587.000 -568.000 1.000 -2.556 588.000 9002.000 ysig <- data@left
N <- length(ysig)
x <- seq(1,N,1)
T <- 6.85
freq <- x/T
length(freq)## [1] 66150 # sample 66150 3 sec rate 22050
# N = len(X) python
# n = np.arange(N)
# T = N/sr
# freq = n/T
# amplitude data left / right speaker
# increments start at 1
time <- T/x
ysig <- ysig / 2^(data@bit -1)
plot(time,ysig,xlim=c(0,0.002))
length(time)## [1] 66150 par(mar=c(1,1,1,1)) # figure margin error fix
ffty <- fftw::FFT(ysig) # discrete fourier transform
amp <- abs(ffty) # DFT is a sequence of complex sinusoids
# find all harmonics with fft()
plot(freq,amp/N,type="l",xlim = c(0,2000))
FT and Signal Averaging
Use R to generate noisy sine waves and use signal averaging to see reduction in S/N ratio.
a useful functionin R is random can be used to add noise to a signal. try with exponential function and added noise. try with large noise that is 10% (SD) of the initial signal.
plot average signal
after 1 after 10 after 100
how does the average noise change with n (=number of signals)