Analysis of Popping Popcorn Audio File

Application of Big Data Techniques

Audio processing workshop:

Analysis of popping corn audioclips

as an exercise in Big Data analysis

CS424 - Big Data

Friday, 8 September 2023

13:30 - 15:00

:::

ANALYSIS OF AUDIO FILES

• Amplitude
• Sampling rate
• Tracks

Audo sampling rate:

plot(-100,0,xlim=c(0,3),ylim=c(0,96000),xlab="",ylab="Frequency",xaxt="n",
     main="Maximun frequency detected")
text(0.75,0,"Preception Limits")
text(2,0,"Device specification")

points(0.75,36000,col="blue",pch=19)
text(0.75,36000,"36,000 Hz",pos=2)
text(0.75,36000,"Teenagers",pos=4)

points(0.75,22000,col="blue",pch=19)
text(0.75,22000,"22,000 Hz",pos=2)
text(0.75,22000,"Adults",pos=4)

points(0.75,14500,col="blue",pch=19)
text(0.75,14500,"14,500 Hz",pos=2)
text(0.75,14500,"Elderly",pos=4)

points(2,8000,col="blue",pch=19)
text(2,8000,"8000 Hz",pos=2)
text(2,8000,"Telephone systems",pos=4)

points(2,16000,col="blue",pch=19)
text(2,16000,"16000 Hz",pos=2)
text(2,16000,"Cassette Tape recordings",pos=4)

points(2,48000,pch=19,col="blue")
text(2,48000,"48000 Hz",pos=2)
text(2,48000,"Music CD",pos=4)

points(2,96000,pch=19,col="blue")
text(2,96000,"96000 Hz",pos=2)
text(2,96000,"Recording studio track",pos=4)

Sound waves

library(tuneR)

Warning: package 'tuneR' was built under R version 4.2.3

zz = sine(440)
x = 1:20000
y = 128*sin(x/19) + 128* sin(x/67) + 128*sin(x/43)
zz@left = y

plot(x/1000,zz@left,type="l",xlab="time in sec",
     ylab="amptitude",main="Sampled 3 Frequencies Modulation") 

Fast Fourier Transform

library(seewave)

Warning: package 'seewave' was built under R version 4.2.3

 spectro(zz,flim=c(0,1),channel=1)

Audio files as a data vector

• Audio recordings can be multitrack
• Typically recorded at 44-1000 KHz

Reading audio files

TuneR library

• Process audio files:
  • generate or read/write audio filessome data: such as sine
  • read or write audio files: (such as readWave, writeWave)
  • represent or construct channelled Wave files: (such as Wave, WaveMC)
  • transform Wave objects (bind, channel, downsample, extractWave, mono, stereo)
  • play Wave objects (play)
• Analysing signals
  • calculate periodograms of a signal (periodogram, Wspec)
  • estimate the fundamental frequencies (FF, FFpure)
• Convert sound to music
  • derive the corresponding notes and melody (noteFromFF, melodyplot)
  • quantization and plot corresponding music (quantize,quantplot)
  • create sheet music lilyinput

library(tuneR)

Loading a sound clip as a data vector

SoundData <- tuneR::readWave('popcorn-popping.wav')
summary(SoundData)

Wave Object
    Number of Samples:      2425500
    Duration (seconds):     55
    Samplingrate (Hertz):   44100
    Channels (Mono/Stereo): Stereo
    PCM (integer format):   TRUE
    Bit (8/16/24/32/64):    24

Summary statistics for channel(s):

          Min. 1st Qu. Median      Mean 3rd Qu.    Max.
left  -4843869  -14542      0  -8.73294   15368 4625660
right -4183111  -17246      0 -24.29296   17544 4820805

play(SoundData)

Visualizing the sound clip

plot(SoundData)

Seewave

• sound wave analysis, manipulation, display, editing and synthesis
• processes oscillograms, spectral content, frequency coherence, dominant frequency, analytic signal, 2D and 3D spectrograms.
• resonance quality factor, entropy, cross correlation and autocorrelation, zero-crossing,
• Homepage: homepage

Analysis of popping corn audio file

• Download the sound file from Freesound

Popping popcorn

Run a spectro analysis

library(seewave)
 spectro(SoundData,flim=c(0,10))

Filter out the relevant frequencies

d = ffilter(SoundData,channel=1,from=0,to=20)
spectro(d,f=44000, channel=1,flim=c(0,0.75))

Check the amplitude of the filtered sound clip

plot(c(1:length(d))/44000,d,type="l",
     main="Filtered sound clip",
     xlab="seconds", ylab="Amplitude")

Analyze the pop of a typical kernel pop

plot(c(100000:130000)/44000,d[100000:130000],type="l",
     main="Filtered sound clip",xlab="seconds")
text(2.44,300, "Kernel cover rips",pos=4,col="red")
rect(2.433,-390,2.445,360,col="#99000066",border="red")
text(2.46,130, "Steam and starch escapes",col="blue",pos=4)
rect(2.45,-100,2.49,100,col="#00009966",border="blue")
text(2.65,-70,"Pot echo/reverbation",col="purple")
rect(2.5,-50,2.8,50,col="#ff00ff66",border="purple")

Rectified wave form

dRec = abs(d[100000:130000])
plot(c(100000:130000)/44000,dRec,type="l",main="Recified sound clip")
text(2.44,300, "Kernel cover rips",pos=4,
     col="red",cex=1.2)
rect(2.431,0,2.443,360,col="#99000066",border="red")
text(2.46,130, "Steam and starch escapes",
     col="blue",pos=4,cex=1.2)
rect(2.45,0,2.49,100,col="#00009966",border="blue")
text(2.65,70,"Pot echo/reverbation",
     col="purple",cex=1.2)
rect(2.5,0,2.8,50,col="#ff00ff66",border="purple")

Smoothed rectified wave form

dSmo = rep(0,1500)
for( i in 0:1499) {
   imn = 20*i
   dSmo[i+1] = max(c(dRec[c(imn:(imn+19)) + 1]))
}
plot(2.273 +c(1:1500)/2200,dSmo,type="l",
      main="Smoothed rectified sound clip",
    xlab="seconds",ylab="Amplitude")
text(2.44,300, "Kernel cover rips",pos=4,
     col="red",cex=1.2)
rect(2.431,0,2.443,360,col="#99000066",border="red")
text(2.46,130, "Steam and starch escapes",
     col="blue",pos=4,cex=1.2)
rect(2.45,0,2.49,100,col="#00009966",border="blue")
text(2.65,70,"Pot echo/reverbation",
     col="purple",cex=1.2)
rect(2.5,0,2.8,50,col="#ff00ff66",border="purple")

Select a threshold for the pop indicator

colors1 =c("skyblue","darkred")
colors2 =c("black","red")

plot(2.273 +c(1:1500)/2200,dSmo,type="l",
    main="Smoothed rectified sound clip",
    xlab="seconds",ylab="Amplitude")
for (i in 2:1500) {
  ci = 1+ (dSmo[i]>60)
     lines(2.273 +c((i-1):i)/2200,dSmo[(i-1):i],
    col=colors2[ci],lwd=2)

  points(2.273 +c(i)/2200,60,
          col=colors1[ci],pch=19,cex=0.5)
  if (dSmo[i]>60) {
    lines(2.273 +c(i,i)/2200,c(60,dSmo[i]),col="red")
  }
}

Create a database of all pops

ddRec = abs(ffilter(d,f=44000,from=0,to=0.8))
ddSmoLen = length(ddRec)/20
ddSmo = rep(0,ddSmoLen)

for(i in 1:ddSmoLen) {
  start = (i-1)*20
  ddSmo[i] = max(ddRec[1+ (start:(start+19))])
}
plot((1:length(ddSmo))/2200,ddSmo,type="l")
lines(c(1,length(ddSmo))/2200,c(.20,.20),col="red",lwd=2)

pkht = 0
pktotal =0
state = 0
peaks = c()
for (i in 1:length(ddSmo)) {
  if ((state == 0) && (ddSmo[i] > .060)) {
    pkht = ddSmo[i]
    pkstart = i
    pktotal = i
    state = 1
  } else if (state == 1) {
    
       if (ddSmo[i] > .060) {
          pktotal = pktotal + ddSmo[i]
          if (ddSmo[i] > pkht) { pkht = ddSmo[i] }
       } else {
          peaks = rbind(peaks,
                  c(pkstart,i-1,i - pkstart, pkht,pktotal))
         # cat(pkstart, i-1, pkstart-i, pkht,"\n")
          state = 0
       }
  }

}

Area

plot(peaks[,4],peaks[,5],pch=19,col="#00009933",
     main="Peak width vs Total Peak Area",
     xlab="Peak width",
     ylab="Area under the peak")

plot(peaks[,3],peaks[,5],col="#00009933",
     main="Max Peak height vs Peak Area",
     xlab="Max peak height",ylab="Peak Area",pch=19)

Set threshold for kernel popping

pops = peaks[peaks[,4]>0.6,]
hist(pops[,5],br=30,main="Frequency distribution of peak area",
     xlab="Peak area",col="blue")

Show popping behavior

pb = hist(pops[,1]/2200,xlab="seconds",
     main="Popping behavior",br=30)

Count the pops

There are 243 popped kernels of corn.

Classify the pops

k = kmeans(pops,3)
colorz=c("red","orange","darkgreen")
plot(pops[,4],pops[,3], col=colorz[k$cluster],pch=19)

Fit a surge function to model the pops per second

\[(pop/sec) = \frac{A \cdot t^C}{e^{Bt}}\] \[ \log \left(\frac{pops}{sec}\right) = \log(A)+ C \log(t) - B t\]

Modelling

pps = log10(pb$counts)
t = pb$mids
logt = log10(pb$mids)
lm = lm(pps ~ logt + t)

pps = log10(pb$counts[5:22])
t = pb$mids[5:22]
logt = log10(pb$mids[5:22])
lm2 = lm(pps ~ logt + t)


pps = log10(pb$counts[7:21])
t = pb$mids[7:21]
logt = log10(pb$mids[7:21])
lm3 = lm(pps ~ logt + t)

plot(pb$mids,pb$counts,col="blue",type="b",ylim=c(0,26))
predicta = 10^(lm$coefficients[1] + lm$coefficients[2]*log10(pb$mids) + lm$coefficients[3] * pb$mids) 
predictb = 10^(lm2$coefficients[1] + lm2$coefficients[2]*log10(pb$mids) + 
                 lm2$coefficients[3] * pb$mids) 
predictc = 10^(lm3$coefficients[1] + lm3$coefficients[2]*log10(pb$mids) + 
                 lm3$coefficients[3] * pb$mids) 
                 
lines(pb$mids, predicta, col="purple")
lines(pb$mids, predictb, col="red")
lines(pb$mids, predictc, col="cyan")
text(pb$mids,pb$counts,1:23,pos=3,col="blue")

legend(10,25,c("A 1-23","B 5-22","C 7-21"),fill=c("purple","red","cyan"))

Determine the mean absolute deviation (MAD) of the prediction.

\[MAD = \sum_{i=1}^n \frac{|\hat y_i - y_i|}{n}\]

• \(MAD_A =\) 4.7673064

• \(MAD_B =\) 3.0490908

• \(MAD_C =\) 2.4443408

Predict the popping cycle of a pot with 50 or 100 kernels of corn given the same conditions.

SoundData <- tuneR::readWave('try1.wav')
summary(SoundData)

Wave Object
    Number of Samples:      4890112
    Duration (seconds):     110.89
    Samplingrate (Hertz):   44100
    Channels (Mono/Stereo): Stereo
    PCM (integer format):   TRUE
    Bit (8/16/24/32/64):    16

Summary statistics for channel(s):

        Min. 1st Qu. Median       Mean 3rd Qu.  Max.
left  -32768    -152      0 -0.6494784     151 29744
right -32768    -149      0 -0.6406131     148 29282

colors= c("blue","red")
rdata = abs(SoundData@left)
plot(rdata,type="l",col="#00009933",lwd=2)
max = length(rdata)
unit = 500
maxlen = max /unit
rd2 = rep(0,maxlen)

for (i in 1: maxlen) {
  rd2[i] = max(rdata[1+ c(((i-1)*unit):(i*unit))])
  
}
lines(c(1:maxlen)*unit,rd2,col="#99000033",lwd=2)

hist(rdata,br=20)

hist(rd2,br=20)

Counting

state = 1
count = 0
pops =c()
popstart= 0
popheight = 0
poparea = 0
for (i in 1:maxlen){
   if (state ==1) {
      if (rd2[i] > 4000) {
        state = 2
        popstart = i
        poparea = poparea + rd2[i]
      } 
     } else if (state == 2) {
        if (rd2[i] > 4000) {
          poparea = poparea + rd2[i]
          if (rd2[i] > popheight) {
            popheight = rd2[i]
          }
        } else {
          state = 1
          pops = rbind(pops,c(count,popstart,popheight,poparea,i-popstart))
          count = count + 1
        }   
     }
}  

Estimated number of pops: 44

Popping behavior

hist(pops[,2]/440,br=15,main="Popping behavior",xlab="Seconds")

Document the experimental results

• Title: A descriptive name for this research
• Abstract: A synopsis of the work and outcome
• Introducti`on: What prompted this work? What do we intend to discover?
• Methodology W:hat research methodology will we use?
• Results: What did we find out?
• Discussion: What do the findings suggest?
• References: What citations and related work are there?

Bibliography

• Sam Carcangno, 2013. Basic sound processing in R
• Dena J. Clink, 2022. Tutorial in R audio processing
• Uwe Ligges et als, 2023. {tuneR}: Analysis of Music and Speech
• J. Sueur, T. Aubin and C. Simonis, 2008. Seewave: a free modular tool for sound analysis and synthesis Bioacoustics, vol 18, pgs 213-226

Possible Class Projects

• Effect of butter and different vegetable oils on popping behavior at a given power setting
• Effect of power settings on popping behavior
• Comparison of the popping behavior of different brands of popcorn
• Correlation of popcorn flavor and texture at different stages of the popping cycle
• Identification of the optimal conditions for popping popcorn