Using MATLAB (779505), страница 51
Текст из файла (страница 51)
For multi-input systems, each row of a matrixcorresponds to a sample point with each input appearing as columns of thematrix.The functiony = filter(b,a,x)processes the data in vector x with the filter described by vectors a and b,creating filtered data y.The filter command can be thought of as an efficient implementation of thedifference equation. The filter structure is the general tapped delay-line filterdescribed by the difference equation below, where n is the index of the currentsample, na is the order of the polynomial described by vector a and nb is theorder of the polynomial described by vector b. The output y(n), is a linearcombination of current and previous inputs, x(n) x(n-1) ..., and previousoutputs, y(n-1) y(n-2) ...a ( 1 )y ( n ) = b ( 1 )x ( n ) + b ( 2 )x ( n – 1 ) + … + b ( nb )x ( n – nb + 1 )– a ( 2 )y ( n – 1 ) – … – a ( na )y ( n – na + 1 )Suppose, for example, we want to smooth our traffic count data with a movingaverage filter to see the average traffic flow over a 4-hour window.
This processis represented by the difference equation1111y ( n ) = --- x ( n ) + --- x ( n – 1 ) + --- x ( n – 2 ) + --- x ( n – 3 )4444The corresponding vectors area = 1;b = [1/4 1/4 1/4 1/4];13-40Difference Equations and FilteringNote Enter the format command, format rat, to display and enter datausing the rational format.Executing the commandload count.datcreates the matrix count in the workspace.For this example, extract the first column of traffic counts and assign it to thevector x.x = count(:,1);The 4-hour moving-average of the data is efficiently calculated withy = filter(b,a,x);Compare the original data and the smoothed data with an overlaid plot of thetwo curves.t = 1:length(x);plot(t,x,'-.',t,y,'-'), grid onlegend('Original Data','Smoothed Data',2)13-4113Data Analysis and Statistics120Original DataSmoothed Data1008060402000510152025The filtered data represented by the solid line is the 4-hour moving average ofthe observed traffic count data represented by the dashed line.For practical filtering applications, the Signal Processing Toolbox includesnumerous functions for designing and analyzing filters.13-42Fourier Analysis and the Fast Fourier Transform (FFT)Fourier Analysis and the Fast Fourier Transform (FFT)Fourier analysis is extremely useful for data analysis, as it breaks down asignal into constituent sinusoids of different frequencies.
For sampled vectordata, Fourier analysis is performed using the discrete Fourier transform(DFT).The fast Fourier transform (FFT) is an efficient algorithm for computing theDFT of a sequence; it is not a separate transform. It is particularly useful inareas such as signal and image processing, where its uses range from filtering,convolution, and frequency analysis to power spectrum estimation.This section:• Summarizes the Fourier transform functions• Introduces Fourier transform analysis with an example about sunspotactivity• Calculates magnitude and phase of transformed data• Discusses the dependence of execution time on length of the transformFunction SummaryMATLAB provides a collection of functions for computing and working withFourier transforms.FFT Function SummaryFunctionDescriptionfftDiscrete Fourier transform.fft2Two-dimensional discrete Fourier transform.fftnN-dimensional discrete Fourier transform.ifftInverse discrete Fourier transform.ifft2Two-dimensional inverse discrete Fourier transform.ifftnN-dimensional inverse discrete Fourier transform.absMagnitude.13-4313Data Analysis and StatisticsFFT Function Summary (Continued)FunctionDescriptionanglePhase angle.unwrapUnwrap phase angle in radians.fftshiftMove zeroth lag to center of spectrum.cplxpairSort numbers into complex conjugate pairs.nextpow2Next higher power of two.IntroductionFor length N input sequence x, the DFT is a length N vector, X.
fft and ifftimplement the relationshipsNX(k)=∑x ( n )en–1– j2π ( k – 1 ) ------------- N 1≤k≤Nn=1N1x ( n ) = ---N∑X ( k )en–1j2π ( k – 1 ) ------------- N 1≤n≤Nk=1Note Since the first element of a MATLAB vector has an index 1, thesummations in the equations above are from 1 to N. These produce identicalresults as traditional Fourier equations with summations from 0 to N-1.If x(n) is real, we can rewrite the above equation in terms of a summation ofsine and cosine functions with real coefficientsN1x ( n ) = ---Nk=113-442π ( k – 1 ) ( n – 1 )2π ( k – 1 ) ( n – 1 )- + b ( k ) sin -------------------------------------------∑ a ( k ) cos -----------------------------------------NNFourier Analysis and the Fast Fourier Transform (FFT)wherea ( k ) = real ( X ( k ) ), b ( k ) = – imag ( X ( k ) ), 1 ≤ n ≤ NFinding an FFTThe FFT of a column vector xx = [4 3 7 -9 1 0 0 0]’ ;is found withy = fft(x)which results iny =6.000011.4853-2.0000-5.485318.0000-5.4853-2.000011.4853- 2.7574i-12.0000i+11.2426i-11.2426i+12.0000i+ 2.7574iNotice that although the sequence x is real, y is complex.
The first componentof the transformed data is the constant contribution and the fifth elementcorresponds to the Nyquist frequency. The last three values of y correspond tonegative frequencies and, for the real sequence x, they are complex conjugatesof three components in the first half of y.Example: Using FFT to Calculate Sunspot PeriodicitySuppose, we want to analyze the variations in sunspot activity over the last 300years. You are probably aware that sunspot activity is cyclical, reaching amaximum about every 11 years. Let’s confirm that.Astronomers have tabulated a quantity called the Wolfer number for almost300 years.
This quantity measures both number and size of sunspots.Load and plot the sunspot data.load sunspot.datyear = sunspot(:,1);13-4513Data Analysis and Statisticswolfer = sunspot(:,2);plot(year,wolfer)title('Sunspot Data')Sunspot Data2001801601401201008060402001700175018001850190019502000Now take the FFT of the sunspot data.Y = fft(wolfer);The result of this transform is the complex vector, Y. The magnitude of Ysquared is called the power and a plot of power versus frequency is a“periodogram.” Remove the first component of Y, which is simply the sum of thedata, and plot the results.N = length(Y);Y(1) = [];power = abs(Y(1:N/2)).^2;nyquist = 1/2;freq = (1:N/2)/(N/2)*nyquist;plot(freq,power), grid onxlabel('cycles/year')title('Periodogram')13-46Fourier Analysis and the Fast Fourier Transform (FFT)Periodogram72x 101.81.61.41.210.80.60.40.2000.050.10.150.20.250.3cycles/year0.350.40.450.5The scale in cycles/year is somewhat inconvenient.
Let’s plot in years/cycle andestimate what one cycle is. For convenience, plot the power versus period(where period = 1./freq) from 0 to 40 years/cycle.period = 1./freq;plot(period,power), axis([0 40 0 2e7]), grid onylabel('Power')xlabel('Period(Years/Cycle)')13-4713Data Analysis and Statistics72x 101.81.61.4Power1.210.80.60.40.200510152025Period(Years/Cycle)30In order to determine the cycle more precisely,[mp,index] = max(power);period(index)ans =11.076913-483540Fourier Analysis and the Fast Fourier Transform (FFT)Magnitude and Phase of Transformed DataImportant information about a transformed sequence includes its magnitudeand phase.
The MATLAB functions abs and angle calculate this information.To try this, create a time vector t, and use this vector to create a sequence xconsisting of two sinusoids at different frequencies.t = 0:1/100:10-1/100;x = sin(2*pi*15*t) + sin(2*pi*40*t);Now use the fft function to compute the DFT of the sequence. The code belowcalculates the magnitude and phase of the transformed sequence. It uses theabs function to obtain the magnitude of the data, the angle function to obtainthe phase information, and unwrap to remove phase jumps greater than pi totheir 2*pi complement.y = fft(x);m = abs(y);p = unwrap(angle(y));Now create a frequency vector for the x-axis and plot the magnitude and phase.f = (0:length(y)-1)'*100/length(y);subplot(2,1,1), plot(f,m),ylabel('Abs.
Magnitude'), grid onsubplot(2,1,2), plot(f,p*180/pi)ylabel('Phase [Degrees]'), grid onxlabel('Frequency [Hertz]')The magnitude plot is perfectly symmetrical about the Nyquist frequency of 50hertz. The useful information in the signal is found in the range 0 to 50 hertz.13-4913Data Analysis and Statistics500Abs. Magnitude40030020010000102030405060708090100102030405060Frequency [Hertz]70809010042.5x 10Phase [Degrees]21.510.50−0.50FFT Length Versus SpeedYou can add a second argument to fft to specify a number of points n for thetransformy = fft(x,n)With this syntax, fft pads x with zeros if it is shorter than n, or truncates it ifit is longer than n.
If you do not specify n, fft defaults to the length of the inputsequence.The execution time for fft depends on the length of the transform. It is fastestfor powers of two. It is almost as fast for lengths that have only small primefactors. It is typically several times slower for lengths that are prime or whichhave large prime factors.The inverse FFT function ifft also accepts a transform length argument.For practical application of the FFT, the Signal Processing Toolbox includesnumerous functions for spectral analysis.13-5014Function FunctionsFunction Summary .
. . . . . . . . . . . . . . . . 14-3Representing Functions in MATLAB . . . . . . . . . 14-4Plotting Mathematical Functions. . . . . . . . . . 14-6Minimizing Functions and Finding Zeros . . . . .Minimizing Functions of One Variable . . . . . . . .Minimizing Functions of Several Variables . . . . . .Setting Minimization Options .
. . . . . . . . . .Finding Zeros of Functions . . . . . . . . . . . . .Tips . . . . . . . . . . . . . . . . . . . . . .Troubleshooting . . . . . . . . . . . . . . . . .Converting Your Optimization Code to MATLAB Version5 Syntax . . . . . . . . . . . . . . . . . ......... 14-9. 14-914-1014-1114-1214-1514-15. 14-16Numerical Integration (Quadrature) . . . . . . . . 14-19Example: Computing the Length of a Curve . . .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
