模擬與數(shù)字信號

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

　　一個連續(xù)的正弦信號。X為時間，Y為每天因為潮汐引起的水位高度

　　數(shù)字化后的信號，每一個點代表采樣點

　　Aliasing

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

　　一個更接近于真實的模擬信號源

　　如果我們對它采樣，大概會得到。。。這樣：

　　用一條線把采樣點連接起來，采集的到的數(shù)字波形可以看出明顯的上下振動。由于高頻噪聲和原來的低頻真實信號相疊加，最后采集出來的數(shù)據(jù)和原來的數(shù)據(jù)相比很難看。這是因為采樣的頻率遠低于噪聲信號的頻率，Aliasing發(fā)生了。

　　避免Aliasing

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

　　Nyquist Theorem

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

　　Anti-Aliasing Filters

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

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

　　RMS

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

　　讓我們手工用matlab擼一個rms函數(shù)：建一個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)生一些隨機數(shù)，然后計算他們的RMS:

rms(rand(1,1000))

　　FFT 傅里葉變換

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

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

　　所謂fft 就是把X軸從時間換成頻率. ifft,就是把X軸坐標從頻率變成時間

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

　　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.

　　然后搞一個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

　　下面準備開始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

　　如果仔細觀察，發(fā)現(xiàn)定點的值是50，而我們sin波形的正弦信號幅度是1啊。。這是matlab fft返回結(jié)果定義的問題，我們只需要除以用于fft運算的采樣點個數(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軸頂點整成1啦。那么X軸怎么辦呢，下面把X軸整成Hz...

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

　　frequency resolution= sample rate / points in FFT

　　頻率分辨率 = 采樣率 / 一次輸入FFT轉(zhuǎ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，美美噠

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

　　%plot first half of data

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

　　xlabel('Frequency (Hz)');

　　ylabel('Amplitude');

　　好了，大功告成

　　測量某一個頻率成分信號的大小

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

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

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

　　[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)濾波器)

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

　　FIR filter

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

　　在這個濾波器中，可以看到每次把前三個數(shù)據(jù)進行平均(分別乘以0.33333)然后輸出。這三個系數(shù)的不同組合(0.3333, 0.333, 0.3333)就組成了各種FIR濾波器。這些系數(shù)叫做filter coefficients. 或者叫做taps. 在matlab中，他們叫做 b. 下面是一個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點moving average 濾波器效果

　　End Bit 問題

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

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

　　總之，我們的5點滑動濾波的計算公式是：

　　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ù)要比原始信號"遲緩" 很多。。這叫做相位延時

　　通過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點moving average filter會使高頻分量衰減。這是一種low pass filter, 通低頻，阻高頻

　　小知識

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

　　ntaps=20;

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

　　IIR filter

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

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

　　先產(chǎn)生一個采樣頻率1000Hz， 2000個采樣點的隨機信號，然后用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有個神奇叫做 濾波器設(shè)計工具，以GUI的方式搞出你想要的濾波器。簡直。。哎。沒啥可說的，這玩意不要教程。。

　　IIR 濾波器原理

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

　　輸入信號為x(m). 對應(yīng)的乘以每一個b系數(shù)。在5點moving average filter的例子中，這個b 就是 0.2,0.2,0.2,0.2,0.2. 輸出為y(m). 這樣每一個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)...

　　自動信號檢測

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

　　能量檢測器

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

　　對信號進行濾波，將有用信號提取出來

　　對信號進行整流

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

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

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

　　% 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');

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

　　step1: 對信號進行濾波, 將有用信號提取出來， 我們搞一個 帶寬為10Hz的 帶通濾波器，中心頻點設(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: 整流信號，其實就是取絕對值，把負的信號弄成正的:

　　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: 對整流過的信號進行包絡(luò)。求整流過的信號的包絡(luò)?？梢酝ㄟ^一個低通濾波器實現(xiàn)。這里我們用一個6極點butter 濾波器實現(xiàn)，截止頻率是 10Hz。過了這個濾波器，濾波后的信號基本就成了原信號的包絡(luò)。注意信號的起始位置，因為過了濾波器，所有可以看出一點點的相位滯后

　　% 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è)置閾值，搞一個變量gated, 當原信號超過一個閾值的時候，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ù)，這個函數(shù)類似于求倒數(shù)，新的值=現(xiàn)在的值-之前一次的值：

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

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

　　圖像處理

　　圖形實際上可以看做三重矩陣，因為每個像素是RGB三個顏色分量，每個像素都用三個值描述，所以是三重矩陣。讓我們先搞一張最簡單的，使用向量建立起來的圖像：

　　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

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

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

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

　　讀取和操作圖像

　　matlab讀取和操作圖像很簡單

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

　　image(r)

　　示例圖片，一張海洋珊瑚的航拍圖, 綠色的是珊瑚，藍色的是海洋

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

　　r(1,1,1:3)

　　ans(:,:,1) =

　　66

　　ans(:,:,2) =

　　153

　　ans(:,:,3) =

　　174

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

　　給圖片加閾值

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

　　step1: 把藍色分量取出來，繪制灰度圖

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

　　%mage(r)

　　rb= r(:,:,3);

　　imagesc(rb);

　　colormap(gray);

　　下面我們對rb進行閾值處理，把<130的標成0，>130的標成1

　　rbt=rb<130; %threshold dark values

　　imagesc(rbt)

　　colormap(gray)

　　最后就簡單了，統(tǒng)計1的個數(shù)，然后除以整個像素數(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就是時域上的坐標。。

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

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

　　FFT輸出的第一個點為：

　　計算FFT的第一個點

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

　　參考

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