Download
1 / 41
410 likes | 422 Views
Gain insights into signals through Fourier analysis in MATLAB. Learn the nuances of time series, sawtooth waveforms, aliasing, spectrum visualization, and DFT. Understand interpolation and analysis of sound signals with practical examples. Explore the challenges of periodic extensions, Gibbs phenomenon, multi-channel data handling, and noise introduction in signals. Discover the importance of FFT and its inverse function for signal processing and computations. Dive into complex forms, discrete time series, and their implications in mathematical analysis and signal processing.
E N D
Math.375 V- Interpolation 3>Sounds and Fourier<>analysis< Vageli Coutsias
A telephone number aka “name that tune”
% 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”.
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
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)
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
%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])
Real form Continuous time forms Complex form
Real form Complex form Discrete time forms
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
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
FFTrecur is of order DFT is of order CSInterp is of order
% 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))
% 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))
% 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))
for i=1:m+1 if (E(i)>.1) i-1 end end x(1:4096); i = 78, 135 1
a = 6000 n = 4096 i = 78,135 2
a = 13000 n = 4096 i = 78, 122 3
a = 19000 n = 4096 i =70, 122 4
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
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.
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
%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);
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
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
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')
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.
[...] = 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
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.
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.
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);
% 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
165,112 a =14000,n =4096 2 65,124 a =41000,n =4096 3 65,137 a =65000,n =4096
4 72,112 a=134000,n=4096 5 72,124 a=134000,n =4096 6 72,137 a =163000,n =4096
7 79,112 a =198000,n =4096 8 79,124 a =231500,n =4096 9 79,137 a =267500,n =4096
* 87,112 a =345000,n =4096 0 87,124 a =299500,n =4096 # 87,137 a =363000,n =4096
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
112 124 137 n=4096 1 2 3 4 5 6 7 8 9 * 0 # 65 72 79 87 The matrix
