CS1371 Introduction to Computing for Engineers

1
CS1371Introduction to Computing for Engineers

• Processing Sound

2
Sound Basics in Matlab
Learning Objectives Learn about how sound is
stored, and how to access and manipulate the data

• Topics
• Anatomy of a sound
• Sound representation in Matlab
• Getting sound into and out of Matlab
• Sound analysis and processing

3
Background

• The physical phenomenon of sound is detected as
  pressure fluctuations on a surface
• ear drum or a microphone.
• Simple sounds are easy to visualize a sine wave
  with certain characteristics would be perceived
  as a pure tone.
• The characteristics of an individual sound
  component are its amplitude and frequency
• In practice, complex sounds are the mathematical
  sum of the contributions of different sound
  components whose amplitude and frequency vary
  continuously, and even discontinuously with time.

4
Sound Amplitude

• Usually reported in decibels
• intensity of a sound in decibels is calculated
  as
• IDB 10 log10( I / I0 )
• where I is the measured pressure fluctuation,
• and I0 is a reference pressure
• usually established as the lowest pressure
  fluctuation a really good ear can detect, 2?10-4
  dynes / cm.2
• a not-too-painful sound level of, say, 100db has
  a pressure amplitude 1010 (10 billion) times as
  loud as the lowest perceivable sound.

5
Sound Frequency

• The normal human ear can detect sounds between
  about 50Hz and 15kHz
• High fidelity systems must be able to record and
  play back sound in these frequency ranges.
• Telephone conversations sound odd because the
  telephone system limits the upper frequency to
  4kHz.
• Excessive exposure to loud noise causes your ears
  to lose sensitivity in the upper and lower
  frequency ranges.

6
Basic digital recording and playback
Sampling frequency
Fs
A / D
D / A
Sound pressure fluctuations
Sound pressure fluctuations
Analog voltage
Sampled signal
Sampled signal
Analog voltage
Digital voltage
Digital voltage
7
Matlab Capabilities

• Matlab has tools for reading sound files in two
  formats
• .wav files (wavread()), and
• .au files (auread()).
• Both return three variables
• vector of sound values in the range (0 1),
• sampling frequency in Hz (samples per second),
  and
• number of bits used to record the data (8 or 16).
• Examine the following example comparing two
  recordings at 16 bit and 8 bit resolution
  respectively.
• gtgt b16, F16, n wavread('brooklyn16.wav')
• gtgt fprintf('bit size of b16 is d\n', n)
• bit size of b16 is 16
• gtgt sound(b16, F16)
• gtgt b8, F8, n wavread('brooklyn8.wav')
• gtgt fprintf('bit size of b8 is d\n', n)
• bit size of b8 is 8
• gtgt sound(b8, F8)

8
Basic digital sound processing
Sampling frequency
Read from data file
Fs
D / A
sound()
Digital Processing
wavread()
9
Types of Processing

• Time domain
• Amplification
• Vector manipulation
• slicing and dicing
• Reversing
• Playback change sampling frequency
• Forcing a change in sample frequency
• Removing data elements
• Frequency Domain
• the Fourier Transform
• Inserting / deleting sounds
• filtering

10
Slicing and Concatenating Sound

• constructing an audio ransom note
• choosing and assembling words from published
  speeches.
• strictly experimental in nature, but
• can produce some interesting results.
• There are a number of web sites that offer free
  .wav files containing memorable speech fragments
  from movies, TV shows and cartoons1.
• 1 beware, however, when looking for materials.
  Some web sites store .wav files in a compressed
  form that Matlab cannot decode. When you find a
  file, save it to disk and make sure that Matlabs
  wavread() function can open it correctly.

11
Slicing and Concatenating Sound

• assembling parts of these speeches into a
  semi-coherent conversation.
• Consider the following scripts
• Step 1 find the problem speech
• houston, Fsh wavread('a13prob.wav')
• sound(houston, Fsh)
• Step 2 extracting the problem speech
• subplot(1, 2, 1)
• plot(houston)
• clip 110000
• prob houston(clipend)2
• subplot(1, 2, 2)
• plot(prob)

Selected experimentally by plotting and listening
Make it louder
12
Slicing and Concatenating Sound

• Complete speech

Last part extracted
13
Extract from Gone with the Wind

• remove my dear from the frankly, my dear ..
  speech, reducing its amplitude by a half
• Step 3 removing my dear
• damn, Fsd wavread('givdamn2.wav')
• subplot(1, 2, 1)
• plot(damn)
• lo 4500
• hi 8700
• sdamn damn(1lo) damn(hiend) .5
• subplot(1, 2, 2)
• plot(sdamn)

Selected experimentally by plotting and listening
Make it softer
14
Slicing and Concatenating Sound

• Complete speech

my dear removed
15
Assembling the Speech

• sound data are stored as column vectors
• concatenation requires semicolons
• Step 4 assembling the speech
• truth, Fst wavread('truth1.wav')
• speech prob sdamn truth .7
• plot(speech)
• sound(speech, Fst)

16
Playing it Backwards

• Since the sound is a vector of data, reversing
  the vector will play it backwards
• Step 5 reversing the speech
• backwards speech(end-11)
• sound( backwards, Fst)

17
Changing Sound Pitch (the hard way)

• manipulate piano.wav to produce a snippet of
  music.
• This file is a recording1 of a single note
  played on a piano.
• If a sound is played at twice its natural
  frequency, it is heard as one musical octave
  higher.
• The steps from one note to the next higher octave
  are divided 12 half note steps.
• These 12 half steps are logarithmically divided
  where the frequency change is 21/12.
• Practical issues
• Cannot concatenate notes because they are played
  at different frequencies
• Duration of the notes must be managed
• 1 credit http//chronos.ece.miami.edu/dasp/sam
  ples/samples.html

18
Code for playing a scale

• note, Fs wavread('piano.wav')
• duration floor(length(note)/8)
• wait duration / Fs
• whole 2(1/6)
• half 2(1/12)
• for i 18
• dur(i) duration
• freq(i) Fs
• sound(note(1duration), Fs)
• pause(wait)
• if (i 3) (i 7)
• Fs Fs half
• duration round(duration half)
• else
• Fs Fs whole
• duration round(duration whole)
• end
• end

19
Code for playing a simple tune

• pitch length
• smallworld 1 3 C
• 1 1 C
• 3 2 E
• 1 2 C
• 2 3 D
• 2 1 D
• 2 4 D
• l length(smallworld)
• for i 1l
• f freq(smallworld(i,1))
• d dur(smallworld(i,1)) smallworld(i,2)
• sound(note(1d), f)
• pause(waitsmallworld(i,2))
• end

20
Changing Sound Frequency (the right way)

• The problem with the previous technique is that
  you are changing the playback frequency.
  Consequently, you cant assemble the whole tune
  into a file that can be saved and played back by
  an ordinary music player.
• We need a technique for changing the frequency of
  a note without changing its playback frequency.
• We can accomplish this by adding or removing data
  samples to extend or shorten the note vector.

half 2(1/12) note Fs wavread(piano.wav)
sound(note, Fs) plays the note at its natural
pitch sound(note(round(1half2end)), Fs)
raises a whole
step sound(note (round(1half-4end)), Fs)
lowers 2 whole steps
21
Script to play a scale

• note, Fs wavread('piano.wav')
• half 2(1/12) whole half2
• steps 0
• for index 18
• sound(note(round(1halfstepsend)), Fs)
• pause(.2)
• if (index 3) (index 7)
• steps steps 1
• else
• steps steps 2
• end
• end

22
Building a Tune

• To make a wav file containing a tune, we have to
  have a description of the song to play, and a
  file with the instrument we wish to use.
• We will describe the song in a (n ? 2) array
  where the first column specifies the note on the
  keyboard, and the second column specifies its
  duration in beats.

key duration
smallworld 1 3 1 1
3 2 1 2 2 3
2 1 2 4
Key 1 2 3 4 5 6 7 8 ½ steps 0 2
4 5 7 9 11 12
23
Building a Tune

• We then write a program that builds the tune by
  putting the notes at the right frequency into the
  right place in a vector

note get
key 1 1 3 get ½ steps
0 0 4 change length
note 1 note 2 note 3 copy to vector
get length 3 1 2 move down
vector beat 0.2 sec ---note 1--- note
2 --note 3--
24
Script to play a tune

• steps 0 2 4 5 7 9 11 12
• col 1 pitch 2 duration
• music 1 3 2 1 3 3 1 1
• 3 2 1 2 3 4 2 3
• 3 1 4 1 4 1 3 1
• 2 1 4 8 3 3 4 1
• 5 3 3 1 5 2 3 2
• 5 4 4 3 5 1 6 1
• 6 1 5 1 4 1 6 4

note, Fs wavread(instrument) half
2(1/12) ln length(music) s
sum(music) totalW s(2) totalSize totalW
beat Fs
length(note) tune zeros(totalSize, 1) start
1 for stp 1ln pitch music(stp, 1)
pow steps(pitch) bit note(round((1(half
pow)end)) lb length(bit) tune(start
(start lb - 1) ) bit start start
beat music(stp, 2) Fs end sound(tune,
Fs) wavwrite(tune, Fs, outName)
25
Front-end for flexibility
function synthesizer( music, beat, outName,
instrument, steps ) if nargin 0 clear
clc warning off all music 1 3 1 1 3
2 1 2 2 3 2 1 2 4 end if nargin lt 5
use the major scale steps 0 2 4 5 7 9 11
12 if nargin lt 4 instrument
'piano.wav' if nargin lt 3
outName 'default.wav' if nargin lt
2 beat 0.2 end
end end end
26
Test Program
clear clc minor 0 2 3 5 7 8 10 12 major 0
2 4 5 7 9 11 12 doremi 1 3 2 1 3 3 1 1 3
2 1 2 3 4 2 3 3 1 4 1 4 1 3 1
2 1 4 8 3 3 4 1 5 3 3 1 5 2 3 2
5 4 4 3 5 1 6 1 6 1 5 1 4 1 6
4 synthesizer all defaults pause(4) synthesize
r(doremi, .25, 'doremiPiano.wav' ) change the
output file name pause(15) synthesizer(doremi,
.15, 'doremiViolin.wav', ...
'violin.wav' ) change the
instrument pause(10) synthesizer(doremi, .2,
'doremiTrumpetMinor.wav', ...
'tpt.wav', minor ) minor key
27
Questions?
28
Fourier Transforms
Dynamic display of sound energy vs frequency in
blocks
Equalizer changes the gain in selected frequency
bands
29
Fourier Transform

• Converts sound vector
• from a sampled time history
• to frequency domain representation

Fourier Transform
Inverse Fourier Transform
30
Fast Fourier Transform

• Common implementation for Fourier Transform
• High efficiency extracting all frequencies
• Specific relationships between time and frequency
  domains

fft()
tmax
?t
ifft()
npts same for each representation
npts same for each representation
?f
fmax
(Not necessarily drawn to scale)
31
FFT Relationships

• Time domain is a vector of npts
• Points are type double in range /- 1 (for
  sound() playback)
• Sampled at ?t 1/Fs
• tmax npts ?t
• Frequency domain is also a vector of npts
• Points are type complex symmetrical complex
  conjugates
• Max frequency at npts/2
• Real parts are mirrored with same sign
• Imaginary parts are mirrored with opposite sign
• Sampled at ?f 1/ tmax
• fmax npts ?f/2

32
Characteristics of a Transform

• create 10 sec of an 8 Hz sine wave, plot the
  first second, do the fft, and examine the
  spectrum.
• subplot(2, 2, 1)
• dt 1/400 sampling period (sec)
• pts 10000 number of points
• f 8 frequency
• t (1pts) dt time array for plotting
• x sin(2pift)
• df 1 / t(end) the frequency interval
• fmax df pts / 2
• plot(t(1end/25), x(1end/25))
• title('time domain sine wave') xlabel('time
  (sec)') ylabel('amplitude')
• subplot(2,2,3)
• Y fft(x) perform the transform
• f (1pts) 2 fmax / pts frequencies for
  plotting
• plot(f, real(Y))
• title('real part') xlabel('frequency (Hz)')
  ylabel('energy')
• subplot(2,2,4)
• plot(f, imag(Y))

33
Resulting Plots
1 sec of a 10 sec recording
Spectrum is mirrored about its center fmax
occurs here really
Real and imaginary components at 8Hz
Real and imaginary components at 8Hz
34
Typical Frequency Domain Processing

• Once in the frequency domain, we can do the
  following things
• Analyze the content of certain sounds, e.g.
  musical instruments
• Add sounds in the frequency domain
• Perform digital filtering

35
Analyzing sounds

• Function to display a music sample from the
  University of Miamis Audio and Signal Processing
  Laboratory.1
• All the instruments are carefully playing a
  steady note at 250 Hz
• function inst(name, ttl) 
• x, Fs wavread(name '.wav')
• N length(x)
• dt 1/Fs sampling period (sec)
• t (1N) dt time array for plotting
• Y abs(fft(x)) perform the transform
• mx max(Y)
• Y Y 100 / mx
• df 1 / t(end) the frequency interval
• fmax df N / 2
• f (1N) 2 fmax / N
• up floor(N/10)
• plot(f(1up), Y(1up) )
• title(ttl) xlabel('frequency (Hz)')
  ylabel('energy')
• 1 http//chronos.ece.miami.edu/dasp/samples/sam
  ples.html

36
Analyzing sounds

• Script to show music samples
• rows 4
• cols 2
• subplot(rows, cols, 1)
• inst('sax','Saxophone')
• subplot(rows, cols, 2)
• inst('flute', 'Flute')
• subplot(rows, cols, 3)
• inst('tbone','Trombone')
• subplot(rows, cols, 4)
• inst('piano', 'Piano')
• subplot(rows, cols, 5)
• inst('tpt', 'Trumpet')
• subplot(rows, cols, 6)
• inst('mutetpt', 'Muted Trumpet')
• subplot(rows, cols, 7)
• inst('violin', 'Violin')
• subplot(rows, cols, 8)
• inst('cello', 'Cello')

37
Resulting plot
38
Add sounds in the frequency domain

• building a signal with two sine waves (3 Hz and 8
  Hz) in the time domain, performing the fft,
  adding a 50Hz wave in the spectrum and using the
  inverse transform (ifft()) to reconstruct and
  play the new signal
• rows 2
• cols 2
• set up two sine waves
• dt 1/400 sampling period (sec)
• N 10000 number of points
• f1 3 first frequency
• f2 8 second frequency
• t (1N) dt time array for plotting
• x sin(2pif1t) sin(2pif2t)
• subplot(rows, cols, 1)
• plot(t(1(1/dt)), x(1(1/dt)))
• title('original signal') xlabel('time (sec)')

39
Add sounds (continued)

• df 1 / t(end) the frequency interval
• fmax df N / 2
• Y fft(x) perform the transform
• f (1N) 2 fmax / N frequencies for
  plotting
• subplot(rows, cols, 2) plot the imaginary
  part
• plot(f(1N/3), imag(Y(1N/3)))
• title('imag spectrum') xlabel('frequency (Hz)')
• find the maximum value and location
• level1, index1 max(imag(Y(1(N/2))).2)
• level2, index2 max(imag(Y((N/21)end)).2)
• newF 50
• newBasic 50 / df
• newI1 newBasic 1
• newI2 N - newBasic 1
• newV complex(0, -sqrt(level1)/2)
• Y(newI1) newV
• Y(newI2) -newV
• subplot(rows, cols, 4) plot the imaginary
  part
• plot(f(1N/3), imag(Y(1N/3)))

40
Resulting plot
41
Perform digital filtering

• In general, digital filtering is complex both
  mathematically and computationally.
• We can contrive a simple exercise by comparing
  the spectra for the trumpet and muted trumpet.
• We should be able to take the trumpet spectrum,
  remove the high-frequency components and
  reproduce the muted trumpet sound.
• This kind of filtering only works well when the
  spectral energy is zero in the areas where
  changes are made.
• If the spectrum is not zero where we insert
  changes, the step functions introduced in the
  frequency domain produce loud clicks or ringing
  in the time domain.

42
Digitally muting a trumpet

• x, Fs wavread('tpt.wav')
• N length(x)
• sound(x, Fs)
• dt 1/Fs sampling period (sec)
• t (1N) dt time array for plotting
• pause(t(N))
• Y fft(x) perform the transform
• df 1 / t(end) the frequency interval
• fmax df N / 2
• f (1N) 2 fmax / N
• subplot(1, 2, 1)
• frac floor(N/10)
• plot(f(1frac), abs(Y(1frac)))
• title('trumpet original')
• xlabel('frequency (Hz)')
• ylabel('energy')

43
Digitally muting a trumpet (continued)

• subplot(1, 2, 2)
• fhalf 925 Hz
• ihalf floor(fhalf / df)
• foff 1100 Hz
• ioff floor(foff / df)
• Y(ihalfend-ihalf) Y(ihalfend-ihalf) 0.2
• Y(ioffend-ioff) 0
• plot(f(1frac), abs(Y(1frac)))
• title('trunc trump') xlabel('freq (Hz)')
  ylabel('energy')
• y abs(ifft(Y))
• mx max(y)
• y y / mx
• disp('synthetic muted')
• sound(y, Fs)
• pause(t(N))
• disp('real muted')
• x, Fs wavread('mutetpt.wav')
• mx max(x)
• x x / mx

44
Muted Trumpet Plot
45
Summary
Action Items Review the lecture Review the
material in our textbook Think of other sound
manipulations to implement and see if you can
create a function to do this.
Learning Matlab can store and manipulate arrays
that represent sound
46
Questions?
47
(No Transcript)