MATLAB课程：代码示例之Mathematics（一）

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

MATLAB课程：代码示例之Mathematics（一）

Using FFT

This example shows how to use the FFT function to analyze the variations in sunspot activity over the last 300 years.

Sunspot activity is cyclical, reaching a maximum about every 11 years. Let's confirm that. Here is a plot of a quantity called the Zurich sunspot relative number, which measures both number and size of sunspots. Astronomers have tabulated this number for almost 300 years.

load sunspot.datyear = sunspot(:,1);relNums = sunspot(:,2);plot(year,relNums)title('Sunspot Data')

Here is a closer look at the first 50 years.

plot(year(1:50),relNums(1:50),'b.-');

The fundamental tool of signal processing is the FFT, or fast Finite Fourier Transform. To take the FFT of the sunspot data type the following.

The first component of Y, Y(1), is simply the sum of the data, and can be removed.

Y = fft(relNums);Y(1) = [];

A graph of the distribution of the Fourier coefficients (given by Y) in the complex plane is pretty, but difficult to interpret. We need a more useful way of examining the data in Y.

plot(Y,'ro')title('Fourier Coefficients in the Complex Plane');xlabel('Real Axis');ylabel('Imaginary Axis');

The complex magnitude squared of Y is called the power, and a plot of power versus frequency is a "periodogram".

n = length(Y);power = abs(Y(1:floor(n/2))).^2;nyquist = 1/2;freq = (1:n/2)/(n/2)*nyquist;plot(freq,power)xlabel('cycles/year')title('Periodogram')

The scale in cycles/year is somewhat inconvenient. We can plot in years/cycle and estimate the length of one cycle.

plot(freq(1:40),power(1:40))xlabel('cycles/year')

Now we plot power versus period for convenience (where period = 1./freq). As expected, there is a very prominent cycle with a length of about 11 years.

period = 1./freq;plot(period,power);axis([0 40 0 2e+7]);ylabel('Power');xlabel('Period (Years/Cycle)');

Finally, we can fix the cycle length a little more precisely by picking out the strongest frequency. The red dot locates this point.

hold on;index = find(power == max(power));mainPeriodStr = num2str(period(index));plot(period(index),power(index),'r.', 'MarkerSize',25);text(period(index)+2,power(index),['Period = ',mainPeriodStr]);hold off;

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：mathematics Mathematic MATLAB课程 Thematic MATLAB MATLAB课程 代码示例 Mathematics MATLAB