关于MATLAB 做FFT的两个问题

%%这两个都是MATLAB中的例子
% 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.dat
year=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;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

例2

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); 
y = x + 2*randn(size(t));     % Sinusoids plus noise
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;为什么要除以L
f = Fs/2*linspace(0,1,NFFT/2+1);这个和f=fs/N*k是一致的吗?

% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')为什么作图时是2*abs,标记的是abs?

 

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值