Math'375 V Interpolation 3 Sounds and Fourier analysis

1
Math.375 V- Interpolation 3gtSounds and
Fourierltgtanalysislt

• Vageli Coutsias

2
A telephone number
aka name that tune
3
A single dial pulse
4
script sounds1 y,Fs,bitswavread('tele_tones_1
.wav') xy(7441000) szlength(x) td
linspace(0,sz,sz)' plot(td,x) sound(x,Fs)
Caution! This is not a periodic signal. We needed
to chop it in order to analyze it as a periodic
time series. But chopping introduces its
own Gibbs phenomenon.
5
f is a periodic function
In practise, a time series of values of f is
available a periodic extension of f is assumed
so that if the f-values at the two end-points are
different, the function behaves as if it had a
jump there. A Gibbs phenomenon results introducing
noisein the computed spectrum
6
For multichannel data, each column of a matrix
represents one channel. Each row of such a matrix
then corresponds to a sample point. A
three-channel signal y that consists of x, 2x,
and x/pi is y x 2x x/pi
Y(1)
Y(2)
Y(3)
7
begin with a vector representing a time base.
Consider generating data with a 1000 Hz
sample frequency. An appropriate time vector
is t (00.0011)' Given t
you can create a sample signal y consisting of
two sinusoids, one at 50 Hz and one at 120 Hz
with twice the amplitude. y sin(2pi50t)
2sin(2pi120t) The new variable y, formed
from vector t, is also 1001 elements long. You
can add normally distributed white noise to the
signal and graph the first fifty points using
8
a sample signal with its time vector t
(00.0011)' y sin(2pi50t)
2sin(2pi120t) randn('state',0) yn y
0.5randn(size(t)) plot(t(150),yn(150))
z
t t.2 square(4t)
a sawtooth wave of prescribed properties fs
10000 t 01/fs1.5 x sawtooth(2pi50t) pl
ot(t,x), axis(0 0.2 -1 1)
9
Real form
Continuous time forms
Complex form
10
Real form
Complex form
Discrete time forms
11
The DFT and its inverse
12
function y DFT(x,isign) isign/-1 1
point sample -gt freq. domain coeffs. -1
freq.dom.coef.-gt point values n length(x) y
x(1)ones(n,1) if ngt1 mn/2
wexp(isignpisqrt(-1)/m) v w.(0n-1)'
for k 2n z rem((k-1)(0n-1)',n )1
y y v(z)x(k) end if isign 1
y y/m end end
13
function y FFTRecur(x) n2k, x col. n
length(x) if n 1 y x else m n/2
yT FFTRecur(x(12n)) yB
FFTRecur(x(22n)) d exp(-2pisqrt(-1)/n).
(0m-1)' z d.yB y yTz yT-z
end
14
FFTrecur is of order
DFT is of order
CSInterp is of order
15
script interface y,Fs,bitswavread('tele_tones
_1.wav') ainput('starting point') ninput('lengt
h of series (power of 2)') xy(aan) mn/2 y
fft(x(1n)) P2 struct('a',real(y(1m1))/m,'b
', -imag(y(2m))/m) P2.a(1)P
2.a(1)/2P2.a(m1)P2.a(m1)/2 Bvals
linspace(0,n/Fs,n1)' Gvals CSeval(P2,n/Fs,Bval
s) plot(Bvals(1100),Gvals(1100)-x(1100))
16
script aliasing z,Fs,bitswavread('tele_tones_
1.wav') ainput('starting point') ninput('length
of series (power of 2)') strideinput('stride
(power of 2) lt length)') xz(astridean)
mn/2/stride y fft(x(1n/stride)) P2
struct('a',real(y(1m1))/m,'b',-imag(y(2m))/m)
P2.a(1)P2.a(1)/2P2.a(m1)P2.a(m1)/2 Bvals
linspace(0,n/Fs,n1)' Gvals CSeval(P2,n/Fs,Bval
s) plot(Bvals(1100stride),Gvals(1100stride)

-z(1100stride))
17
script spectrum y,Fs,bitswavread('tele_tones_
1.wav') ainput('starting point') ninput('length
of series (power of 2)') xy(aan) mn/2)
sound(x,Fs) y fft(x(1n)) P2
struct('a',real(y(1m1))/m,'b',
-imag(y(2m))/m) P2.a(1)P2.a(1)/2P2.a(m1)P
2.a(m1)/2 E(2m)sqrt(P2.a(2m).2P2.b(1m-1).
2) E(1)abs(P2.a(1))E(m1)abs(P2.a(m1)) plot(
E(1500))
18
for i1m1 if (E(i)gt.1) i-1
end end
x(14096) i 78, 135
1
19
a 6000 n 4096 i 78,135
2
20
a 13000 n 4096 i 78, 122
3
21
a 19000 n 4096 i 70, 122
4
22
a 25500 n 4096 i 78, 135
a 38000 n 4096 i 78, 135
5
7
XXY-ZXXX
a 32000 n 4096 I 78,135
6
23
APPENDIX a
FFT Discrete Fourier transform. FFT(X) is the
discrete Fourier transform (DFT) of vector X.
For matrices, the FFT operation is applied to
each column. For N-D arrays, the FFT
operation operates on the first non-singleton
dimension. FFT(X,N) is the N-point FFT,
padded with zeros if X has less than N points
and truncated if it has more. FFT(X,,DIM)
or FFT(X,N,DIM) applies the FFT operation
across the dimension DIM. See also
IFFT, FFT2, IFFT2, FFTSHIFT.
24
For length N input vector x, the DFT is a length
N vector X, with elements
The inverse DFT
(computed by IFFT) is given by

25
an example of use FFT, DFT, FFTRecur n512mn/2
t (0(1/n)1)' kklinspace(0,512,513)' y
sin(2pi50t) 2cos(2pi120t) z1fft(y(1n))
-----------------------FFT x1real(ifft(z1)) z
2FFTRecur(y(1n)) --------------FFTRecur P2str
uct('a',real(z2(1m1))/m,'b',imag(z2(2m))/m) P2
.a(1)P2.a(1)/2P2.a(m1)P2.a(m1)/2 x2
CSeval(P2,1,t) z3DFT(y(1n),1)
------------------DFT P3struct('a',real(z3(1m1
))/m,'b',imag(z3(2m))/m) P3.a(1)P3.a(1)/2P3.a(
m1)P3.a(m1)/2 x3 CSeval(P3,1,t)
26
APPENDIX b
HELP Reference\MatlabFunctionReference\ File I/O
\Audio and Audio/Video General
SPARCstation-Specific Sound Functions
Microsoft WAVE Sound Functions Audio Video
Interleaved (AVI) Functions Microsoft Excel
Functions Lotus123 Functions
27
Microsoft WAVE Sound
Functions wavplay Play sound on PC-based
audio output device wavread Read Microsoft
WAVE (.wav) sound file wavrecord Record
sound using PC-based audio input device wavwrite
Write Microsoft WAVE (.wav) sound file
28
Wavread Read Microsoft WAVE (.wav) sound
file Graphical Interface As an alternative to
auread, use the Import Wizard. To activate the
Import Wizard, select Import Data from the File
menu. Syntax y
wavread('filename') y,Fs,bits
wavread('filename') ...
wavread('filename',N) ...
wavread('filename',N1 N2) ...
wavread('filename','size')
29
Description wavread supports multichannel data,
with up to 16 bits per sample. y
wavread('filename') loads a WAVE file specified
by the string filename, returning the sampled
data in y. The .wav extension is appended if no
extension is given. Amplitude values are in the
range -1,1. y,Fs,bits wavread('filename')
returns the sample rate (Fs) in Hertz and the
number of bits per sample (bits) used to encode
the data in the file. ... wavread('filename',
N) returns only the first N samples from each
channel in the file.
30
... wavread('filename',N1 N2) returns
only samples N1 through N2 from each channel
in the file. siz wavread('filename','size')
returns the size of the audio data contained in
the file in place of the actual audio data,
returning the vector siz samples
channels. See Also auread, wavwrite
31
Wavwrite Write Microsoft WAVE (.wav) sound
file Syntax wavwrite(y,'filename')
wavwrite(y,Fs,'filename')
wavwrite(y,Fs,N,'filename') Description wavwrite
supports multi-channel 8- or 16-bit WAVE
data. wavwrite(y,'filename') writes a WAVE file
specified by the string filename. The data
should be arranged with one channel per column.
Amplitude values outside the range -1,1
are clipped prior to writing. Wavwrite(y,Fs,'filen
ame') specifies the sample rate Fs, in
Hertz, of the data. Wavwrite(y,Fs,N,'filename')
forces an N-bit file format to be written,
where N lt 16.
32
Sound Convert vector into sound Syntax
sound(y,Fs) sound(y)
sound(y,Fs,bits) Description sound(y,Fs) sends
the signal in vector y (with sample frequency
Fs) to the speaker on PC and most UNIX
platforms. Values in y are assumed to be in
the range . Values outside that range are
clipped. Stereo sound is played on platforms
that support it when y is an n-by-2 matrix.
sound(y) plays the sound at the default sample
rate or 8192Hz. sound(y,Fs,bits) plays the
sound using bits number of bits/sample, if
possible. Most platforms support bits 8 or bits
16.
33
Appendix c
p10.m - polynomials and corresponding
equipotential curves (N.Trefethen, Spectral
Methods in Matlab) N 16 clf for i 12
if i1, s 'equispaced points'
x -1 2(0N)/N end if i2, s
'Chebyshev points' x
cos(pi(0N)/N) end p poly(x)
34
Plot p(x) over -1,1 xx -1.0051 pp
polyval(p,xx) subplot(2,2,2i-1)
plot(x,0x,'.','markersize',13), hold on
plot(xx,pp), grid on set(gca,'xtick',-1.51),
title(s) Plot equipotential curves
subplot(2,2,2i) plot(real(x),imag(x),'.','mar
kersize',13), hold on axis(-1.4 1.4 -1.12
1.12) xgrid -1.4.021.4 ygrid
-1.12.021.12 xx,yy meshgrid(xgrid,ygrid
) zz xx1iyy pp polyval(p,zz) levels
10.(-40) contour(xx,yy,abs(pp),levels),
title(s), colormap(1e-61 1 1) end
35
(No Transcript)
36
1 65,112 a 14000,n 4096
2 65,124 a 41000,n 4096
3 65,137 a 65000,n 4096
37
4 72,112 a134000,n4096
5 72,124 a134000,n 4096
6 72,137 a 163000,n 4096
38
7 79,112 a 198000,n 4096
8 79,124 a 231500,n 4096
9 79,137 a 267500,n 4096
39
87,112 a 345000,n 4096
0 87,124 a 299500,n 4096
87,137 a 363000,n 4096
40
frequencies in Hertz
1218.75 1343.75 1468.75
687.5 781.25 843.75 937.5
1 2 3 4 5
6 7 8 9 0
The matrix
41
112 124 137
n4096
1 2 3 4 5
6 7 8 9 0
65 72 79 87
The matrix