模拟与数字信号

　　我们本身生活在一个模拟量的世界里，所谓模拟量，即连续变化量，屋里的温度是连续变化的，时间是连续变化的，诸如此类。而计算机是数字系统，他不能处理模拟量，而只能处理离散量，这意味着我们要把连续的模拟量进行采样，得到一系列离散的数字量。

　　一个连续的正弦信号。X为时间，Y为每天因为潮汐引起的水位高度

　　数字化后的信号，每一个点代表采样点

　　Aliasing

　　所有自然界的信号都不是很干净的，都会有噪声。下面是一个更接近真实的潮汐水位高度随时间变化的数据集：

　　一个更接近于真实的模拟信号源

　　如果我们对它采样，大概会得到。。。这样：

　　用一条线把采样点连接起来，采集的到的数字波形可以看出明显的上下振动。由于高频噪声和原来的低频真实信号相叠加，最后采集出来的数据和原来的数据相比很难看。这是因为采样的频率远低于噪声信号的频率，Aliasing发生了。

　　避免Aliasing

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

　　Nyquist Theorem

　　这个定理告诉我们，如果想要完整采集某个频率的信号，那么你至少要用2倍于该信号的最高频率的采样率来采集，否则 Aliasing就会发生。举个例子：如果你知道潮汐变化最短以一天为周期，那么你至少半天就需要采集一个潮汐水位变化。实际应用中，通常用更高（>2倍）的采样率去采集信号

　　Anti-Aliasing Filters

　　除了提高采样率，另外一种避免Aliasing的方法是使用 Anti-Aliasing Filter. 通常在数字化之前使用一个low-pass filter把噪声滤掉即可。举例：采样之前安装一个低通滤波器，截止频率为10Hz. 那么你只需要一个20Hz的采样率就可以把你感兴趣的信号采集进来。高频噪声在采样之前就被模拟低通滤波器干掉了。

　　但是：一定要在数字化(采样)之前进行低通滤波(模拟低通滤波)。否则如果采样率不够，则必然发生Aliasing, 噪声会混进你感兴趣的低频信号中，这时候再采用数字低通滤波就没啥效果了。当然，如果你采样频率够高，那么采样进来后才进行数字低通滤波也是OK的。而且绝大多数应用就是这么干的。

　　RMS

　　RMS=root,mean,square. 用来描述信号质量。计算方法: 一组数据，先平方，再求均值，最后开根。

　　让我们手工用matlab撸一个rms函数：建一个rms.m 写入以下内容：

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

使用也很简单，产生一些随机数，然后计算他们的RMS:

rms(rand(1,1000))

　　FFT 傅里叶变换

　　傅里叶变换是信号处理里强大的工具，他可以帮助我们分析信号的频率成分，本文不会讲解傅里叶变化的数学原理，而会侧重matlab实践。

　　fft和ifft(逆傅里叶变换)可以把信号进行时域和频域转换如下图。

　　所谓fft 就是把X轴从时间换成频率. ifft,就是把X轴坐标从频率变成时间

　　先用matlab产生一个100Hz正弦信号, 先产生时间序列

　　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函数就可以了，注意fft返回的一组复数，包含了频率成分和相位成分，我们要绝对值一把fft的结果：

　　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

　　如果仔细观察，发现定点的值是50，而我们sin波形的正弦信号幅度是1啊。。这是matlab fft返回结果定义的问题，我们只需要除以用于fft运算的采样点个数的一半即可：

　　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转换的采样点数

　　在我们的例子中， 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的信号怎么右面还有一个对称的？。。这个你就别问了，蜜汁数学。通常我们只取前面一半就好：

　　%plot first half of data

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

　　xlabel('Frequency (Hz)');

　　ylabel('Amplitude');

　　好了，大功告成

　　测量某一个频率成分信号的大小

　　之前说过，FFT变换的频率分辨率取决于信号采样频率和进入FFT函数的数据量。如果你有一个1000点的数据，原信号是8000Hz，那么频率分辨率就是8Hz

　　获得最大赋值对应的频率

　　通过fft很容易让我们知道赋值最大的信号频点在哪里：

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

　　Returns:

　　magnitude = 1

　　index = 3(从0开始，所以2Hz是3)

　　数字滤波器

　　数字滤波器用于一些不需要的频率信号。比如工频噪声(50Hz)正弦信号等。通常有两类滤波器:

　　FIR(有限冲击响应滤波器)

　　IIR(无限冲击响应滤波器)

　　虽然名字听起来很高大上，其实他们计算并不算复杂。本文侧重matlab应用，并不会涉及深奥的理论知识。

　　FIR filter

　　其实你可能之前用过FIR filter只不过没有意识到而已， moving average(滑动平均)滤波器就是FIR滤波器的一种。在一些应用中，有一个窗口，每一次新的数据进来都在窗口进行一次平局然后输出。

　　在这个滤波器中，可以看到每次把前三个数据进行平均(分别乘以0.33333)然后输出。这三个系数的不同组合(0.3333, 0.333, 0.3333)就组成了各种FIR滤波器。这些系数叫做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 问题

　　你可能已经注意到了， 对于moving average filter 一开始滤波的时候没有之前的数据做平均， 这可咋整？matlab的 filter函数是这么整的：

　　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).

　　可以想到，如果滤波器的阶数很大(N) 那么滤波后的数据要比原始信号"迟缓" 很多。。这叫做相位延时

　　通过FFT看滤波后的频率响应

　　我们把原始信号和滤波后信号用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, 通低频，阻高频

　　小知识

　　其实搞一个滑动平均滤波器的系数可以用下面的简便方法：

　　ntaps=20;

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

　　IIR filter

　　IIR和FIR类似，不过进入了反馈机制。即下一次的滤波输出不仅仅和上几次的输入信号有关，还和上几次的输出信号有关。IIR比FIR"效率"更高，通常用更少的系数就可以达到很好的滤波结果，但是IIR也有缺点，由于引入了反馈机制，一些特定的系数组成的IIR滤波器可能不稳定，造成输出结果崩溃。。

　　根据IIR滤波器的不同系数 有很多经典的IIR滤波器，什么Butterworth, Chebyshew, Bessel之类的，反正都是以牛人的名字命名的。本文只讨论一种滤波器：butterworth. 下面还是上例子：

　　先产生一个采样频率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 函数是matlab内置函数，输入截止频率和阶数，返回滤波器系数，b,a。b 就跟FIR 一样的b(前馈系数)， a指的是后馈系数。

　　可以看出同样是5阶滤波，butter filter的滤波效果单从幅频响应来说 比moving average 好了不少。

　　FDAtool

　　matlab有个神奇叫做 滤波器设计工具，以GUI的方式搞出你想要的滤波器。简直。。哎。没啥可说的，这玩意不要教程。。

　　IIR 滤波器原理

　　使用matlab的帮助文档，可以看到一个框图，介绍滤波器是如何工作的：

　　输入信号为x(m). 对应的乘以每一个b系数。在5点moving average filter的例子中，这个b 就是 0.2,0.2,0.2,0.2,0.2. 输出为y(m). 这样每一个x(m)乘上对应的系数然后再加在一起 组成了 y.　　IIR 滤波器引入了'a' 系数反馈环节。每一次滤波，上一次的输出也要程序对应的系数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)...

　　自动信号检测

　　信号检测指的是在带有噪声的信号中发掘有用信号。本文还是侧重于实践，尽量避免理论数学知识

　　能量检测器

　　一般的能量检测器包含以下步骤

　　对信号进行滤波，将有用信号提取出来

　　对信号进行整流

　　对整流过的信号进行包络

　　设置阈值，检测有用信号

　　还是一个例子, 这次我们把一个正弦信号，藏 在噪声信号中的一个随机位置:

　　% 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的 带通滤波器，中心频点设置在需要提取的正弦信号位置，这一波操作代码如下：

　　% 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: 对整流过的信号进行包络。求整流过的信号的包络。可以通过一个低通滤波器实现。这里我们用一个6极点butter 滤波器实现，截止频率是 10Hz。过了这个滤波器，滤波后的信号基本就成了原信号的包络。注意信号的起始位置，因为过了滤波器，所有可以看出一点点的相位滞后

　　% 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: 给信号设置阈值，搞一个变量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函数，这个函数类似于求倒数，新的值=现在的值-之前一次的值：

　　这样我们就能轻松的得到信号开始的时间了：使用 find(d_gated==1)

　　思考题：使用 find函数来找到正弦信号什么时候结束

　　图像处理

　　图形实际上可以看做三重矩阵，因为每个像素是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

　　这个程序创建了一个全是0,和1的矩阵，双击变量y可以看到y是由0和1的列组成的

　　把Y以图像的形式显示出来：

　　第一幅创建的图形

　　读取和操作图像

　　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 。

　　给图片加阈值

　　我们的目的是统计图片里珊瑚所占整个图像的比例。基本思路就是把图像里不是很蓝的部分找出来，标记成1，蓝的的地方找出来，标记成0.然后 用1的个数除以整个图像像素数即可得珊瑚所占比例

　　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)

　　最后就简单了，统计1的个数，然后除以整个像素数即可：

　　% 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的数量

　　例: N=4, 则 n和m为0到3

　　FFT输出的第一个点为：

　　计算FFT的第一个点

　　最后，复数去绝对值的公式

　　参考

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