在處理海洋和大氣時(shí)間序列時(shí)，我們經(jīng)常需要用傅里葉變換來(lái)對(duì)其做頻譜分析，尋找其顯著周期。傅里葉變換的有關(guān)知識(shí)大家可以自行百度，這里我舉一個(gè)有趣的小例子來(lái)說(shuō)明傅里葉變換可能存在的問(wèn)題。下面的代碼改編自MATLAB函數(shù)FFT自帶的例子，有兩個(gè)正弦信號(hào)，一個(gè)頻率是50Hz, 另一個(gè)是150Hz。將這個(gè)兩個(gè)正弦信號(hào)相加，然后做傅里葉變換，結(jié)果如圖1所示，可以看到在50和150Hz處有峰。但是如果這兩個(gè)正弦信號(hào)不是相加，而是相乘呢？結(jié)果如圖2所示，在100和200Hz處有峰。原因很簡(jiǎn)單，大家可以回想起高中學(xué)過(guò)的三角函數(shù)積化和差公式，兩個(gè)正弦信號(hào)相乘可以轉(zhuǎn)化為兩個(gè)正弦信號(hào)相加。傅里葉變換是無(wú)法反映乘法運(yùn)算的，只能反映加法運(yùn)算。在現(xiàn)實(shí)世界里，乘法運(yùn)算是普遍存在的，比如潮波傳播到近海，由于水深變淺，在摩擦項(xiàng)的非線(xiàn)性作用下，會(huì)誕生淺水分潮。

Fs = 1000;??????????????????? % Sampling frequency

T = 1/Fs;???????????????????? % Sample time

L = 1000;???????????????????? % Length of signal

t = (0:L-1)*T;??????????????? % Time vector

x1=sin(2*pi*50*t);

x2=2*sin(2*pi*150*t);

y= x1+x2; %y= x1.*x2;

NFFT = 2^nextpow2(L); % Next power of 2 from length of y

Y = fft(y,NFFT)/L;

f = Fs/2*linspace(0,1,NFFT/2+1);

figure()

plot(f,2*abs(Y(1:NFFT/2+1)))

xlabel('頻率 (Hz)')

ylabel('振幅')

set(gca,'Fontsize',13,'linewidth',1.1)

?

剛才的例子非常簡(jiǎn)單易懂，這里再舉一個(gè)更復(fù)雜的例子。保留頻率為50Hz的正弦信號(hào)，將頻率為150Hz的正弦信號(hào)替換成隨機(jī)噪聲。如果將50Hz的正弦信號(hào)與隨機(jī)噪聲相加，頻譜分析結(jié)果如圖3所示，可以看到在50Hz出有一個(gè)明顯的峰，其他頻率上也有峰，但是都非常?。ù砹穗S機(jī)噪聲）。如果將50Hz的正弦信號(hào)與隨機(jī)噪聲相乘，頻譜分析結(jié)果如圖4所示，在50Hz處已經(jīng)沒(méi)有峰了。我們還是可以用積化和差來(lái)解釋?zhuān)ò央S機(jī)噪聲看成是一系列余弦信號(hào)的疊加）。

Fs = 1000;??????????????????? % Sampling frequency

T = 1/Fs;???????????????????? % Sample time

L = 1000;???????????????????? % Length of signal

t = (0:L-1)*T;??????????????? % Time vector

x1=sin(2*pi*50*t);

x2=0.5*randn(size(t));

y= x1+x2;%y= x1.*x2

NFFT = 2^nextpow2(L); % Next power of 2 from length of y

Y = fft(y,NFFT)/L;

f = Fs/2*linspace(0,1,NFFT/2+1);

figure()

plot(f,2*abs(Y(1:NFFT/2+1)))

xlabel('頻率 (Hz)')

ylabel('振幅')

set(gca,'Fontsize',13,'linewidth',1.1)