模擬與數(shù)字信號

　　我們本身生活在一個(gè)模擬量的世界里,，所謂模擬量,，即連續(xù)變化量，屋里的溫度是連續(xù)變化的,，時(shí)間是連續(xù)變化的,，諸如此類。而計(jì)算機(jī)是數(shù)字系統(tǒng),，他不能處理模擬量,，而只能處理離散量，這意味著我們要把連續(xù)的模擬量進(jìn)行采樣,，得到一系列離散的數(shù)字量,。

　　一個(gè)連續(xù)的正弦信號。X為時(shí)間,，Y為每天因?yàn)槌毕鸬乃桓叨?/p>

　　數(shù)字化后的信號,，每一個(gè)點(diǎn)代表采樣點(diǎn)

　　Aliasing

　　所有自然界的信號都不是很干凈的，都會(huì)有噪聲,。下面是一個(gè)更接近真實(shí)的潮汐水位高度隨時(shí)間變化的數(shù)據(jù)集：

　　一個(gè)更接近于真實(shí)的模擬信號源

　　如果我們對它采樣,，大概會(huì)得到。,。,。這樣：

　　用一條線把采樣點(diǎn)連接起來,，采集的到的數(shù)字波形可以看出明顯的上下振動(dòng)。由于高頻噪聲和原來的低頻真實(shí)信號相疊加,，最后采集出來的數(shù)據(jù)和原來的數(shù)據(jù)相比很難看,。這是因?yàn)椴蓸拥念l率遠(yuǎn)低于噪聲信號的頻率，Aliasing發(fā)生了,。

　　避免Aliasing

　　Aliasing 通常是有害但又不可能100%避免的：

　　Nyquist Theorem

　　這個(gè)定理告訴我們,，如果想要完整采集某個(gè)頻率的信號，那么你至少要用2倍于該信號的最高頻率的采樣率來采集,，否則 Aliasing就會(huì)發(fā)生,。舉個(gè)例子：如果你知道潮汐變化最短以一天為周期，那么你至少半天就需要采集一個(gè)潮汐水位變化,。實(shí)際應(yīng)用中,，通常用更高（>2倍）的采樣率去采集信號

　　Anti-Aliasing Filters

　　除了提高采樣率，另外一種避免Aliasing的方法是使用 Anti-Aliasing Filter. 通常在數(shù)字化之前使用一個(gè)low-pass filter把噪聲濾掉即可,。舉例：采樣之前安裝一個(gè)低通濾波器,，截止頻率為10Hz. 那么你只需要一個(gè)20Hz的采樣率就可以把你感興趣的信號采集進(jìn)來。高頻噪聲在采樣之前就被模擬低通濾波器干掉了,。

　　但是：一定要在數(shù)字化(采樣)之前進(jìn)行低通濾波(模擬低通濾波),。否則如果采樣率不夠，則必然發(fā)生Aliasing, 噪聲會(huì)混進(jìn)你感興趣的低頻信號中,，這時(shí)候再采用數(shù)字低通濾波就沒啥效果了,。當(dāng)然，如果你采樣頻率夠高,，那么采樣進(jìn)來后才進(jìn)行數(shù)字低通濾波也是OK的,。而且絕大多數(shù)應(yīng)用就是這么干的。

　　RMS

　　RMS=root,mean,square. 用來描述信號質(zhì)量,。計(jì)算方法: 一組數(shù)據(jù),，先平方，再求均值,，最后開根,。

　　讓我們手工用matlab擼一個(gè)rms函數(shù)：建一個(gè)rms.m 寫入以下內(nèi)容：

　　function h=rms(data)% h=rms(data)

　　% calculates root-mean-square

　　% of input matrix data

　　% Square the data using cell-by-cell multiplication

　　datasquared=data.*data;

　　% Calculate the mean of the squared data

　　mean_ds=mean(datasquared);

　　% Calculate the square-root of the mean

　　% h is what will be returned to the

　　% calling function

　　h=sqrt(mean_ds);

使用也很簡單，產(chǎn)生一些隨機(jī)數(shù),，然后計(jì)算他們的RMS:

rms(rand(1,1000))

　　FFT 傅里葉變換

　　傅里葉變換是信號處理里強(qiáng)大的工具,，他可以幫助我們分析信號的頻率成分，本文不會(huì)講解傅里葉變化的數(shù)學(xué)原理,，而會(huì)側(cè)重matlab實(shí)踐,。

　　fft和ifft(逆傅里葉變換)可以把信號進(jìn)行時(shí)域和頻域轉(zhuǎn)換如下圖。

　　所謂fft 就是把X軸從時(shí)間換成頻率. ifft,就是把X軸坐標(biāo)從頻率變成時(shí)間

　　先用matlab產(chǎn)生一個(gè)100Hz正弦信號, 先產(chǎn)生時(shí)間序列

　　srate = 100; % 100 Hz sample rate

　　speriod = 1/srate; % sample period (s)

　　dur=1; % duration in seconds

　　t=[0:speriod:dur-speriod];% Create a signal from 0 to 1 s. The sample period is 0.01 seconds.

　　然后搞一個(gè)2Hz的sin波形, 并打印一下

　　freq=2; % frequency of sine wave in Hz

　　s=sin(2*pi*freq.*t); % calculate sine wave.

　　plot(t,s); % plot using t as the time base

　　下面準(zhǔn)備開始fft了,， 使用matlab自帶的fft函數(shù)就可以了,，注意fft返回的一組復(fù)數(shù),，包含了頻率成分和相位成分，我們要絕對值一把fft的結(jié)果：

　　sfft=fft(s); % calculate Fast Fourier Transform

　　sfft_mag=abs(sfft); % take abs to return the magnitude of the frequency components

　　plot(sfft_mag); % plot fft result

　　如果仔細(xì)觀察,，發(fā)現(xiàn)定點(diǎn)的值是50，而我們sin波形的正弦信號幅度是1啊,。,。這是matlab fft返回結(jié)果定義的問題，我們只需要除以用于fft運(yùn)算的采樣點(diǎn)個(gè)數(shù)的一半即可：

　　fftpts=length(s);

　　hpts=fftpts/2; % half of the number points in FFT

　　sfft_mag_scaled=sfft_mag/hpts; % scale FFT result by half of the number of points used in the FFT

　　plot(sfft_mag_scaled); % plot scaled fft result

　　Y軸整好了,。,。

　　好啦，Y軸頂點(diǎn)整成1啦,。那么X軸怎么辦呢,，下面把X軸整成Hz...

　　首先需要知道fft的頻率分辨率為:

　　frequency resolution= sample rate / points in FFT

　　頻率分辨率 = 采樣率 / 一次輸入FFT轉(zhuǎn)換的采樣點(diǎn)數(shù)

　　在我們的例子中， binwidth(頻率分辨率) 就是1Hz.

　　binwidth = srate/fftpts; % 100 Hz sample rate/ 100 points in FFT

　　f=[0:binwidth:srate-binwidth]; % frequency scale goes from 0 to sample rate.Since we are counting from zero, we subtract off one binwidth to get the correct number of points

　　plot(f,sfft_mag_scaled); % plot fft with frequency axis

　　xlabel('Frequency (Hz)'); % label x-axis

　　ylabel('Amplitude'); % label y-axis

　　X,Y軸都搞定了,，趕緊表上X,Y label,，美美噠

　　當(dāng)然我們還有最大的問題沒有解決，明明是2Hz的信號怎么右面還有一個(gè)對稱的,？,。。這個(gè)你就別問了,，蜜汁數(shù)學(xué),。通常我們只取前面一半就好：

　　%plot first half of data

　　plot(f(1:hpts),sfft_mag_scaled(1:hpts));

　　xlabel('Frequency (Hz)');

　　ylabel('Amplitude');

　　好了，大功告成

　　測量某一個(gè)頻率成分信號的大小

　　之前說過,，F(xiàn)FT變換的頻率分辨率取決于信號采樣頻率和進(jìn)入FFT函數(shù)的數(shù)據(jù)量,。如果你有一個(gè)1000點(diǎn)的數(shù)據(jù)，原信號是8000Hz,，那么頻率分辨率就是8Hz

　　獲得最大賦值對應(yīng)的頻率

　　通過fft很容易讓我們知道賦值最大的信號頻點(diǎn)在哪里：

　　[magnitude, index]=max(sfft_mag_scaled(1:hpts))

　　Returns:

　　magnitude = 1

　　index = 3(從0開始,，所以2Hz是3)

　　數(shù)字濾波器

　　數(shù)字濾波器用于一些不需要的頻率信號。比如工頻噪聲(50Hz)正弦信號等,。通常有兩類濾波器:

　　FIR(有限沖擊響應(yīng)濾波器)

　　IIR(無限沖擊響應(yīng)濾波器)

　　雖然名字聽起來很高大上,，其實(shí)他們計(jì)算并不算復(fù)雜。本文側(cè)重matlab應(yīng)用,，并不會(huì)涉及深?yuàn)W的理論知識(shí),。

　　FIR filter

　　其實(shí)你可能之前用過FIR filter只不過沒有意識(shí)到而已， moving average(滑動(dòng)平均)濾波器就是FIR濾波器的一種,。在一些應(yīng)用中,，有一個(gè)窗口，每一次新的數(shù)據(jù)進(jìn)來都在窗口進(jìn)行一次平局然后輸出,。

　　在這個(gè)濾波器中,，可以看到每次把前三個(gè)數(shù)據(jù)進(jìn)行平均(分別乘以0.33333)然后輸出,。這三個(gè)系數(shù)的不同組合(0.3333, 0.333, 0.3333)就組成了各種FIR濾波器。這些系數(shù)叫做filter coefficients. 或者叫做taps. 在matlab中,，他們叫做 b. 下面是一個(gè)moving average filter的 例子:

　　npts=1000;

　　b=[.2 .2 .2 .2 .2]; % create filter coefficients for 5- point moving average

　　x=(rand(npts,1)*2)-1; % raw data from -1 to +1

　　filtered_data=filter(b,1,x); % filter using 'b' coefficients

　　subplot(2,1,1); % 1st subplot

　　plot(x); % plot raw data

　　title('Raw Data');

　　subplot(2,1,2); % 2nd subplot

　　plot(filtered_data); %plot filtered data

　　title('Filtered Data');

　　xlabel('Time')

　　5點(diǎn)moving average 濾波器效果

　　End Bit 問題

　　你可能已經(jīng)注意到了,， 對于moving average filter 一開始濾波的時(shí)候沒有之前的數(shù)據(jù)做平均， 這可咋整,？matlab的 filter函數(shù)是這么整的：

　　filtered_data(2) = x(1)*0.2 + x(2)*0.2.

　　總之,，我們的5點(diǎn)滑動(dòng)濾波的計(jì)算公式是：

　　filtered_data(n) = b(1)*x(n) + b(2)*x(n-1) + b(3)*x(n-2) + b(4)*x(n-3) + b(5)*x(n-4).

　　可以想到，如果濾波器的階數(shù)很大(N) 那么濾波后的數(shù)據(jù)要比原始信號"遲緩" 很多,。,。這叫做相位延時(shí)

　　通過FFT看濾波后的頻率響應(yīng)

　　我們把原始信號和濾波后信號用FFT搞一把:

　　% Perform FFT on original and filtered signal

　　fftpts=npts; % number points in FFT

　　hpts=fftpts/2; % half number of FFT points

　　x_fft=abs(fft(x))/hpts; %scaled FFT of original signal

　　filtered_fft=abs(fft(filtered_data))/hpts; %scaled FFT of filtered signal

　　subplot(2,1,1) %1st subplot

　　plot(x_fft(1:hpts)); %plot first half of data points

　　title('Raw Data');

　　subplot(2,1,2) %2nd subplot

　　plot(filtered_fft(1:hpts));%plot first half data points

　　title('Filtered Data');

　　xlabel('Frequency');

　　可以看出，5點(diǎn)moving average filter會(huì)使高頻分量衰減,。這是一種low pass filter, 通低頻,，阻高頻

　　小知識(shí)

　　其實(shí)搞一個(gè)滑動(dòng)平均濾波器的系數(shù)可以用下面的簡便方法：

　　ntaps=20;

　　b=ones(1,ntaps)/ntaps; % create filter coefficients for ntap-point moving average

　　IIR filter

　　IIR和FIR類似，不過進(jìn)入了反饋機(jī)制,。即下一次的濾波輸出不僅僅和上幾次的輸入信號有關(guān),，還和上幾次的輸出信號有關(guān)。IIR比FIR"效率"更高,，通常用更少的系數(shù)就可以達(dá)到很好的濾波結(jié)果,，但是IIR也有缺點(diǎn)，由于引入了反饋機(jī)制,，一些特定的系數(shù)組成的IIR濾波器可能不穩(wěn)定,，造成輸出結(jié)果崩潰。,。

　　根據(jù)IIR濾波器的不同系數(shù) 有很多經(jīng)典的IIR濾波器,，什么Butterworth, Chebyshew, Bessel之類的，反正都是以牛人的名字命名的,。本文只討論一種濾波器：butterworth. 下面還是上例子：

　　先產(chǎn)生一個(gè)采樣頻率1000Hz,， 2000個(gè)采樣點(diǎn)的隨機(jī)信號，然后用butter 濾一波：

　　% Butterworth IIR Filter

　　srate=1000; % sample rate

　　npts=2000; %npts in signal

　　Nyquist=srate/2; %Nyquist frequency

　　lpf=300; %low-pass frequency

　　order=5; %filter order

　　t=[0:npts-1]/srate; %time scale for plot

　　x=(rand(npts,1)*2)-1; % raw data from -1 to +1

　　[b,a]=butter(order,lpf/Nyquist); %create filter coefficients

　　filtered_data=filter(b,a,x); % filter using 'b' and 'a' coefficients

　　% Calculate FFT

　　fftpts=npts; % number points in FFT

　　hpts=fftpts/2; % half number of FFT points

　　binwidth=srate/fftpts;

　　f=[0:binwidth:srate-binwidth];

　　x_fft=abs(fft(x))/hpts; %scaled FFT of original signal

　　filtered_fft=abs(fft(filtered_data))/hpts; %scaled FFT of filtered signal

　　subplot(2,2,1)

　　plot(t,x);

　　title('Raw Time Series');

　　subplot(2,2,3)

　　plot(t,filtered_data);

　　title('Filtered Time Series');

　　xlabel('Time (s)');

　　subplot(2,2,2)

　　plot(f(1:hpts),x_fft(1:hpts));

　　title('Raw FFT');

　　subplot(2,2,4)

　　plot(f(1:hpts),filtered_fft(1:hpts));

　　title('Filtered FFT');

　　xlabel('Frequency (Hz)');

　　簡單解釋一下程序中可能看不懂的地方：butter 函數(shù)是matlab內(nèi)置函數(shù),，輸入截止頻率和階數(shù),，返回濾波器系數(shù)，b,a,。b 就跟FIR 一樣的b(前饋系數(shù)),， a指的是后饋系數(shù)。

　　可以看出同樣是5階濾波,，butter filter的濾波效果單從幅頻響應(yīng)來說 比moving average 好了不少,。

　　FDAtool

　　matlab有個(gè)神奇叫做 濾波器設(shè)計(jì)工具，以GUI的方式搞出你想要的濾波器,。簡直,。,。哎。沒啥可說的,，這玩意不要教程,。。

　　IIR 濾波器原理

　　使用matlab的幫助文檔,，可以看到一個(gè)框圖,，介紹濾波器是如何工作的：

　　輸入信號為x(m). 對應(yīng)的乘以每一個(gè)b系數(shù)。在5點(diǎn)moving average filter的例子中,，這個(gè)b 就是 0.2,0.2,0.2,0.2,0.2. 輸出為y(m). 這樣每一個(gè)x(m)乘上對應(yīng)的系數(shù)然后再加在一起 組成了 y.　　IIR 濾波器引入了'a' 系數(shù)反饋環(huán)節(jié)。每一次濾波,，上一次的輸出也要程序?qū)?yīng)的系數(shù)a 然后減到本次輸出中：

　　y(n)=b(1)*x(n)+b(2)*x(n-1)+ ... + b(n)*x(n) - a(2)y(n-1) - a(3)*y(n-2)...

　　自動(dòng)信號檢測

　　信號檢測指的是在帶有噪聲的信號中發(fā)掘有用信號,。本文還是側(cè)重于實(shí)踐，盡量避免理論數(shù)學(xué)知識(shí)

　　能量檢測器

　　一般的能量檢測器包含以下步驟

　　對信號進(jìn)行濾波,，將有用信號提取出來

　　對信號進(jìn)行整流

　　對整流過的信號進(jìn)行包絡(luò)

　　設(shè)置閾值,，檢測有用信號

　　還是一個(gè)例子, 這次我們把一個(gè)正弦信號，藏 在噪聲信號中的一個(gè)隨機(jī)位置:

　　% Create Noise Signal and embed sine wave in it at a random location

　　npts=5000; % # points in signal

　　srate=1000; % sample rate (Hz)

　　dur=npts/srate; % signal duration in sec

　　amp_n=3; % noise amplitude

　　amp_t=1; % sine wave amplitude

　　freq=100; % sine wave frequency

　　sinepts=400; % # points in sine wave

　　t=[0:npts-1]/srate; % signal time base

　　sine_t=[0:sinepts-1]/srate; % sine wave time base

　　noise=(rand(1,npts)-.5)*amp_n; %noise signal

　　sinewave=amp_t*(sin(2*pi*freq*sine_t)); % sine wave

　　% random location of sinewave

　　sinepos=floor(rand(1,1)*(npts-sinepts))+1;

　　signal = noise; % make signal equal to noise

　　endsine = sinepos + sinepts-1; %calc index end of sine wave

　　% add sinewave to signal starting at index sinepos

　　signal(sinepos:endsine) = signal(sinepos:endsine) + sinewave;

　　% Plot signal

　　subplot(5,1,1)

　　plot(t,signal);

　　hold on

　　%plot red dot on location of sine wave start

　　plot(sinepos/srate,0,'.r');

　　hold off

　　title('Original Signal');

　　這個(gè)圖把所有每一步的波形全部顯示出來

　　step1: 對信號進(jìn)行濾波, 將有用信號提取出來,， 我們搞一個(gè) 帶寬為10Hz的 帶通濾波器,，中心頻點(diǎn)設(shè)置在需要提取的正弦信號位置，這一波操作代碼如下：

　　% Step 1: Filter signal

　　hbandwidth=5; %half bandwidth

　　nyquist=srate/2;

　　ntaps=200;

　　hpf=(freq-hbandwidth)/nyquist;

　　lpf=(freq+hbandwidth)/nyquist;

　　[b,a]=fir1(ntaps, [hpf lpf]); %calc filter coefficients

　　signal_f=filter(b,a,signal); %filter signal

　　subplot(5,1,2)

　　plot(t,signal_f); %plot filtered signal

　　hold on

　　%plot red dot on location of sine wave

　　plot(sinepos/srate,0,'.r');

　　hold off

　　title('Step 1. Filter Signal');

　　step2: 整流信號,，其實(shí)就是取絕對值,，把負(fù)的信號弄成正的:

　　signal_r=abs(signal_f); %abs value of filtered signal

　　subplot(5,1,3)

　　plot(t,signal_r); %plot filtered signal

　　hold on

　　%plot red dot on location of sine wave

　　plot(sinepos/srate,0,'.r');

　　hold off

　　title('Step 2. Rectify Signal');

　　step3: 對整流過的信號進(jìn)行包絡(luò)。求整流過的信號的包絡(luò),?？梢酝ㄟ^一個(gè)低通濾波器實(shí)現(xiàn)。這里我們用一個(gè)6極點(diǎn)butter 濾波器實(shí)現(xiàn),，截止頻率是 10Hz,。過了這個(gè)濾波器，濾波后的信號基本就成了原信號的包絡(luò),。注意信號的起始位置,，因?yàn)檫^了濾波器，所有可以看出一點(diǎn)點(diǎn)的相位滯后

　　% Step 3: Low-pass filter rectified signal

　　lpf=10/nyquist; % low-pass filter corner frequency

　　npoles=6;

　　[b,a]=butter(npoles, lpf); % Butterworth filter coeff

　　signal_env=filter(b,a,signal_r); %Filter signal

　　subplot(5,1,4)

　　plot(t,signal_env); %plot filtered signal

　　hold on

　　%plot red dot on location of sine wave

　　plot(sinepos/srate,0,'.r');

　　hold off

　　title('Step 3. Envelope Signal');

　　step4: 給信號設(shè)置閾值,，搞一個(gè)變量gated, 當(dāng)原信號超過一個(gè)閾值的時(shí)候,，gated=1, 否則=0

　　% Step 4: Threshold signal

　　threshold=amp_t/2;

　　gated=(signal_env>threshold);

　　subplot(5,1,5)

　　plot(t,signal_env); %plot filtered signal

　　hold on

　　%plot red dot on location of sine wave

　　plot(sinepos/srate,0,'.r');

　　plot(t, gated, 'r'); % plot threshold signal

　　hold off

　　title('Step 4. Threshold Signal');

　　xlabel('Time (s)');

　　setp5：diff一把，找到信號

　　d_gated=diff(gated);

　　plot(d_gated);

　　快使用diff函數(shù),，這個(gè)函數(shù)類似于求倒數(shù),，新的值=現(xiàn)在的值-之前一次的值：

　　這樣我們就能輕松的得到信號開始的時(shí)間了：使用 find(d_gated==1)

　　思考題：使用 find函數(shù)來找到正弦信號什么時(shí)候結(jié)束

　　圖像處理

　　圖形實(shí)際上可以看做三重矩陣，因?yàn)槊總€(gè)像素是RGB三個(gè)顏色分量,，每個(gè)像素都用三個(gè)值描述,，所以是三重矩陣,。讓我們先搞一張最簡單的，使用向量建立起來的圖像：

　　rows=20; % variable for # rows

　　stripes=5;% variable for # stripes

　　% make one column of 1's and an adjacent column of 0's

　　x=horzcat(ones(rows,1),zeros(rows,1));

　　y=[]; %initialize and clear y matrix

　　for n=1:stripes

　　y=horzcat(y,x);% concatenate x onto y

　　end

　　clim=[0 1]; % color limits for imagesc

　　imagesc(y, clim); % plot scaled image

　　colormap(gray); % use gray scale color map

　　這個(gè)程序創(chuàng)建了一個(gè)全是0,和1的矩陣,，雙擊變量y可以看到y(tǒng)是由0和1的列組成的

　　把Y以圖像的形式顯示出來：

　　第一幅創(chuàng)建的圖形

　　讀取和操作圖像

　　matlab讀取和操作圖像很簡單

　　r=imread('reef.jpg','jpg');

　　image(r)

　　示例圖片,，一張海洋珊瑚的航拍圖, 綠色的是珊瑚，藍(lán)色的是海洋

　　通過觀察r可以看出,，r是一個(gè) 三維矩陣,，每一個(gè)維度代表red，green, blue, 比如我們想看看第一個(gè)像素的RGB值,，可以用:

　　r(1,1,1:3)

　　ans(:,:,1) =

　　66

　　ans(:,:,2) =

　　153

　　ans(:,:,3) =

　　174

　　每一個(gè)像素都是由RGB三個(gè)值所組成的,，第一個(gè)像素 red = 66, green = 153, blue=175 。

　　給圖片加閾值

　　我們的目的是統(tǒng)計(jì)圖片里珊瑚所占整個(gè)圖像的比例,?；舅悸肪褪前褕D像里不是很藍(lán)的部分找出來，標(biāo)記成1,，藍(lán)的的地方找出來,，標(biāo)記成0.然后 用1的個(gè)數(shù)除以整個(gè)圖像像素?cái)?shù)即可得珊瑚所占比例

　　step1: 把藍(lán)色分量取出來，繪制灰度圖

　　r=imread('reef.jpg','jpg');

　　%mage(r)

　　rb= r(:,:,3);

　　imagesc(rb);

　　colormap(gray);

　　下面我們對rb進(jìn)行閾值處理,，把<130的標(biāo)成0,，>130的標(biāo)成1

　　rbt=rb<130; %threshold dark values

　　imagesc(rbt)

　　colormap(gray)

　　最后就簡單了，統(tǒng)計(jì)1的個(gè)數(shù),，然后除以整個(gè)像素?cái)?shù)即可：

　　% calculate sum of all white points (=1)

　　reef=sum(sum(rbt)); %reef pts

　　totalpts=prod(size(rbt)); %#pts in image

　　percentreef=reef/totalpts

　　附錄 - DFT(離散傅里葉變換)

　　離散傅里葉變換公式

　　其中

　　X(m) 是變換的輸出X(0), X(1),X(2)...X(N-1). m是輸出頻域上的頻率值. m = 0,1,2,3,4...N-1

　　x(n)是輸入信號x(0), x(1), x(2)... n就是時(shí)域上的坐標(biāo),。。

　　N就是輸入 sample的數(shù)量

　　例: N=4, 則 n和m為0到3

　　FFT輸出的第一個(gè)點(diǎn)為：

　　計(jì)算FFT的第一個(gè)點(diǎn)

　　最后,，復(fù)數(shù)去絕對值的公式

　　參考

　　MATLAB SIMPLIFIED Practical Programming and Signal Processing for Scientists Copyright 2013 David Mann, Loggerhead Instruments