Math.375 V- Interpolation 3 >Sounds and Fourier< >analysis<

1. Math.375 V- Interpolation 3>Sounds and Fourier<>analysis< Vageli Coutsias

2. A telephone number aka “name that tune”

3. A single dial pulse

4. % script sounds1 [y,Fs,bits]=wavread('tele_tones_1.wav'); x=y(744:1000); sz=length(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 “noise”in 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 2*x x/pi] Y(1) Y(2) Y(3)

7. begin with a vector representing atime base. Consider generating data with a 1000 Hz sample frequency. An appropriate time vector is t = (0:0.001:1)'; 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(2*pi*50*t) + 2*sin(2*pi*120*t); 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 = (0:0.001:1)'; y = sin(2*pi*50*t) + 2*sin(2*pi*120*t); randn('state',0); yn = y + 0.5*randn(size(t)); plot(t(1:50),yn(1:50)) % z = [t t.^2 square(4*t)]; % a sawtooth wave of prescribed properties fs = 10000; t = 0:1/fs:1.5; x = sawtooth(2*pi*50*t); plot(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 -> freq. domain coeffs. % -1 freq.dom.coef.-> point values n = length(x); y = x(1)*ones(n,1); if n>1 m=n/2; w=exp(isign*pi*sqrt(-1)/m); v = w.^(0:n-1)'; for k = 2:n z = rem((k-1)*(0:n-1)',n )+1; y = y + v(z)*x(k); end if isign == 1; y = y/m; end end

13. function y = FFTRecur(x) %n=2^k, x col. n = length(x); if n == 1 y = x; else m = n/2; yT = FFTRecur(x(1:2:n)); yB = FFTRecur(x(2:2:n)); d = exp(-2*pi*sqrt(-1)/n).^(0:m-1)'; z = d.*yB; y = [ yT+z ; yT-z ]; end

14. FFTrecur is of order DFT is of order CSInterp is of order

15. % script interface [y,Fs,bits]=wavread('tele_tones_1.wav'); a=input('starting point') n=input('length of series (=power of 2)') x=y(a:a+n); m=n/2; y = fft(x(1:n)); P2 = struct('a',real(y(1:m+1))/m,'b', -imag(y(2:m))/m); P2.a(1)=P2.a(1)/2;P2.a(m+1)=P2.a(m+1)/2; Bvals = linspace(0,n/Fs,n+1)'; Gvals = CSeval(P2,n/Fs,Bvals); plot(Bvals(1:100),Gvals(1:100)-x(1:100))

16. % script aliasing [z,Fs,bits]=wavread('tele_tones_1.wav'); a=input('starting point') n=input('length of series (=power of 2)') stride=input('stride (=power of 2) < length)') x=z(a:stride:a+n); m=n/2/stride; y = fft(x(1:n/stride)); P2 = struct('a',real(y(1:m+1))/m,'b',-imag(y(2:m))/m); P2.a(1)=P2.a(1)/2;P2.a(m+1)=P2.a(m+1)/2; Bvals = linspace(0,n/Fs,n+1)'; Gvals = CSeval(P2,n/Fs,Bvals); plot(Bvals(1:100*stride),Gvals(1:100*stride) -z(1:100*stride))

17. % script spectrum [y,Fs,bits]=wavread('tele_tones_1.wav'); a=input('starting point') n=input('length of series (=power of 2)') x=y(a:a+n); m=n/2;) % sound(x,Fs) y = fft(x(1:n)); P2 = struct('a',real(y(1:m+1))/m,'b', -imag(y(2:m))/m); P2.a(1)=P2.a(1)/2;P2.a(m+1)=P2.a(m+1)/2; E(2:m)=sqrt(P2.a(2:m).^2+P2.b(1:m-1).^2); E(1)=abs(P2.a(1));E(m+1)=abs(P2.a(m+1)); plot(E(1:500))

18. for i=1:m+1 if (E(i)>.1) i-1 end end x(1:4096); 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 n=512;m=n/2; t = (0:(1/n):1)'; kk=linspace(0,512,513)'; y = sin(2*pi*50*t) + 2*cos(2*pi*120*t); z1=fft(y(1:n)); %-----------------------FFT x1=real(ifft(z1)); z2=FFTRecur(y(1:n)); %--------------FFTRecur P2=struct('a',real(z2(1:m+1))/m,'b',imag(z2(2:m))/m); P2.a(1)=P2.a(1)/2;P2.a(m+1)=P2.a(m+1)/2; x2 = CSeval(P2,1,t); z3=DFT(y(1:n),1); %------------------DFT P3=struct('a',real(z3(1:m+1))/m,'b',imag(z3(2:m))/m); P3.a(1)=P3.a(1)/2;P3.a(m+1)=P3.a(m+1)/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 WAVESound Functions Audio Video Interleaved (AVI) Functions Microsoft Excel Functions Lotus123 Functions

27. MicrosoftWAVE 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. WavreadRead 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 MicrosoftWAVE (.wav) sound file Syntaxwavwrite(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,'filename')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 <= 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 = 1:2 if i==1, s = 'equispaced points'; x = -1 + 2*(0:N)/N; end if i==2, s = 'Chebyshev points'; x = cos(pi*(0:N)/N); end p = poly(x);

34. % Plot p(x) over [-1,1]: xx = -1:.005:1; pp = polyval(p,xx); subplot(2,2,2*i-1) plot(x,0*x,'.','markersize',13), hold on plot(xx,pp), grid on set(gca,'xtick',-1:.5:1), title(s) % Plot equipotential curves: subplot(2,2,2*i) plot(real(x),imag(x),'.','markersize',13), hold on axis([-1.4 1.4 -1.12 1.12]) xgrid = -1.4:.02:1.4; ygrid = -1.12:.02:1.12; [xx,yy] = meshgrid(xgrid,ygrid); zz = xx+1i*yy; pp = polyval(p,zz); levels = 10.^(-4:0); contour(xx,yy,abs(pp),levels), title(s), colormap(1e-6*[1 1 1]); end

36. 165,112 a =14000,n =4096 2 65,124 a =41000,n =4096 3 65,137 a =65000,n =4096

37. 4 72,112 a=134000,n=4096 5 72,124 a=134000,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 n=4096 1 2 3 4 5 6 7 8 9 * 0 # 65 72 79 87 The matrix