模擬與數(shù)字信號(hào)
我們本身生活在一個(gè)模擬量的世界里，所謂模擬量，即連續(xù)變化量，屋里的溫度是連續(xù)變化的，時(shí)間是連續(xù)變化的，諸如此類。而計(jì)算機(jī)是數(shù)字系統(tǒng)，他不能處理模擬量，而只能處理離散量，這意味著我們要把連續(xù)的模擬量進(jìn)行采樣，得到一系列離散的數(shù)字量。
一個(gè)連續(xù)的正弦信號(hào)。X為時(shí)間，Y為每天因?yàn)槌毕鸬乃桓叨?/p>
數(shù)字化后的信號(hào)，每一個(gè)點(diǎn)代表采樣點(diǎn)
Aliasing
所有自然界的信號(hào)都不是很干凈的，都會(huì)有噪聲。下面是一個(gè)更接近真實(shí)的潮汐水位高度隨時(shí)間變化的數(shù)據(jù)集：
一個(gè)更接近于真實(shí)的模擬信號(hào)源
如果我們對(duì)它采樣，大概會(huì)得到。。。這樣：
用一條線把采樣點(diǎn)連接起來，采集的到的數(shù)字波形可以看出明顯的上下振動(dòng)。由于高頻噪聲和原來的低頻真實(shí)信號(hào)相疊加，最后采集出來的數(shù)據(jù)和原來的數(shù)據(jù)相比很難看。這是因?yàn)椴蓸拥念l率遠(yuǎn)低于噪聲信號(hào)的頻率，Aliasing發(fā)生了。
避免Aliasing
Aliasing 通常是有害但又不可能100%避免的：
Nyquist Theorem
這個(gè)定理告訴我們，如果想要完整采集某個(gè)頻率的信號(hào)，那么你至少要用2倍于該信號(hào)的最高頻率的采樣率來采集，否則 Aliasing就會(huì)發(fā)生。舉個(gè)例子：如果你知道潮汐變化最短以一天為周期，那么你至少半天就需要采集一個(gè)潮汐水位變化。實(shí)際應(yīng)用中，通常用更高（>2倍）的采樣率去采集信號(hào)
Anti-Aliasing Filters
除了提高采樣率，另外一種避免Aliasing的方法是使用 Anti-Aliasing Filter. 通常在數(shù)字化之前使用一個(gè)low-pass filter把噪聲濾掉即可。舉例：采樣之前安裝一個(gè)低通濾波器，截止頻率為10Hz. 那么你只需要一個(gè)20Hz的采樣率就可以把你感興趣的信號(hào)采集進(jìn)來。高頻噪聲在采樣之前就被模擬低通濾波器干掉了。
但是：一定要在數(shù)字化(采樣)之前進(jìn)行低通濾波(模擬低通濾波)。否則如果采樣率不夠，則必然發(fā)生Aliasing, 噪聲會(huì)混進(jìn)你感興趣的低頻信號(hào)中，這時(shí)候再采用數(shù)字低通濾波就沒啥效果了。當(dāng)然，如果你采樣頻率夠高，那么采樣進(jìn)來后才進(jìn)行數(shù)字低通濾波也是OK的。而且絕大多數(shù)應(yīng)用就是這么干的。
RMS
RMS=root,mean,square. 用來描述信號(hào)質(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);
使用也很簡(jiǎn)單，產(chǎn)生一些隨機(jī)數(shù)，然后計(jì)算他們的RMS:
rms(rand(1,1000))
FFT 傅里葉變換
傅里葉變換是信號(hào)處理里強(qiáng)大的工具，他可以幫助我們分析信號(hào)的頻率成分，本文不會(huì)講解傅里葉變化的數(shù)學(xué)原理，而會(huì)側(cè)重matlab實(shí)踐。
fft和ifft(逆傅里葉變換)可以把信號(hào)進(jìn)行時(shí)域和頻域轉(zhuǎn)換如下圖。
所謂fft 就是把X軸從時(shí)間換成頻率. ifft,就是把X軸坐標(biāo)從頻率變成時(shí)間
先用matlab產(chǎn)生一個(gè)100Hz正弦信號(hào), 先產(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ù)，包含了頻率成分和相位成分，我們要絕對(duì)值一把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波形的正弦信號(hào)幅度是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的信號(hào)怎么右面還有一個(gè)對(duì)稱的？。。這個(gè)你就別問了，蜜汁數(shù)學(xué)。通常我們只取前面一半就好：
%plot first half of data
plot(f(1:hpts),sfft_mag_scaled(1:hpts));
xlabel('Frequency (Hz)');
ylabel('Amplitude');
好了，大功告成
測(cè)量某一個(gè)頻率成分信號(hào)的大小
之前說過，F(xiàn)FT變換的頻率分辨率取決于信號(hào)采樣頻率和進(jìn)入FFT函數(shù)的數(shù)據(jù)量。如果你有一個(gè)1000點(diǎn)的數(shù)據(jù)，原信號(hào)是8000Hz，那么頻率分辨率就是8Hz
獲得最大賦值對(duì)應(yīng)的頻率
通過fft很容易讓我們知道賦值最大的信號(hào)頻點(diǎn)在哪里：
[magnitude, index]=max(sfft_mag_scaled(1:hpts))
Returns:
magnitude = 1
index = 3(從0開始，所以2Hz是3)
數(shù)字濾波器用于一些不需要的頻率信號(hào)。比如工頻噪聲(50Hz)正弦信號(hào)等。通常有兩類濾波器:
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)注意到了， 對(duì)于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ù)要比原始信號(hào)"遲緩" 很多。。這叫做相位延時(shí)
通過FFT看濾波后的頻率響應(yīng)
我們把原始信號(hào)和濾波后信號(hào)用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ù)可以用下面的簡(jiǎn)便方法：
ntaps=20;
b=ones(1,ntaps)/ntaps; % create filter coefficients for ntap-point moving average
IIR filter
IIR和FIR類似，不過進(jìn)入了反饋機(jī)制。即下一次的濾波輸出不僅僅和上幾次的輸入信號(hào)有關(guān)，還和上幾次的輸出信號(hào)有關(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ī)信號(hào)，然后用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)');
簡(jiǎn)單解釋一下程序中可能看不懂的地方：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的方式搞出你想要的濾波器。簡(jiǎn)直。。哎。沒啥可說的，這玩意不要教程。。
IIR 濾波器原理
使用matlab的幫助文檔，可以看到一個(gè)框圖，介紹濾波器是如何工作的：
輸入信號(hào)為x(m). 對(duì)應(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)乘上對(duì)應(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)信號(hào)檢測(cè)
信號(hào)檢測(cè)指的是在帶有噪聲的信號(hào)中發(fā)掘有用信號(hào)。本文還是側(cè)重于實(shí)踐，盡量避免理論數(shù)學(xué)知識(shí)
能量檢測(cè)器
一般的能量檢測(cè)器包含以下步驟
對(duì)信號(hào)進(jìn)行濾波，將有用信號(hào)提取出來
對(duì)信號(hào)進(jìn)行整流
對(duì)整流過的信號(hào)進(jìn)行包絡(luò)
設(shè)置閾值，檢測(cè)有用信號(hào)
還是一個(gè)例子, 這次我們把一個(gè)正弦信號(hào)，藏 在噪聲信號(hào)中的一個(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: 對(duì)信號(hào)進(jìn)行濾波, 將有用信號(hào)提取出來， 我們搞一個(gè) 帶寬為10Hz的 帶通濾波器，中心頻點(diǎn)設(shè)置在需要提取的正弦信號(hào)位置，這一波操作代碼如下：
% 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: 整流信號(hào)，其實(shí)就是取絕對(duì)值，把負(fù)的信號(hào)弄成正的:
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: 對(duì)整流過的信號(hào)進(jìn)行包絡(luò)。求整流過的信號(hào)的包絡(luò)?？梢酝ㄟ^一個(gè)低通濾波器實(shí)現(xiàn)。這里我們用一個(gè)6極點(diǎn)butter 濾波器實(shí)現(xiàn)，截止頻率是 10Hz。過了這個(gè)濾波器，濾波后的信號(hào)基本就成了原信號(hào)的包絡(luò)。注意信號(hào)的起始位置，因?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: 給信號(hào)設(shè)置閾值，搞一個(gè)變量gated, 當(dāng)原信號(hào)超過一個(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一把，找到信號(hào)
d_gated=diff(gated);
plot(d_gated);
快使用diff函數(shù)，這個(gè)函數(shù)類似于求倒數(shù)，新的值=現(xiàn)在的值-之前一次的值：
這樣我們就能輕松的得到信號(hào)開始的時(shí)間了：使用 find(d_gated==1)
思考題：使用 find函數(shù)來找到正弦信號(hào)什么時(shí)候結(jié)束
圖像處理
圖形實(shí)際上可以看做三重矩陣，因?yàn)槊總€(gè)像素是RGB三個(gè)顏色分量，每個(gè)像素都用三個(gè)值描述，所以是三重矩陣。讓我們先搞一張最簡(jiǎn)單的，使用向量建立起來的圖像：
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讀取和操作圖像很簡(jiǎn)單
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);
下面我們對(duì)rb進(jìn)行閾值處理，把<130的標(biāo)成0，>130的標(biāo)成1
rbt=rb<130; %threshold dark values
imagesc(rbt)
colormap(gray)
最后就簡(jiǎn)單了，統(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)是輸入信號(hào)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ù)去絕對(duì)值的公式
參考
MATLAB SIMPLIFIED Practical Programming and Signal Processing for Scientists Copyright 2013 David Mann, Loggerhead Instruments