matlab 数字滤波入门

matlab 数字滤波入门

序言

这是一篇最简单不过的matlab数字信号处理的介绍,里面涉及数字滤波,简单的图像处理和信号检测

1. 时间序列分析入门

模拟与数字信号

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

一个连续的正弦信号。X为时间,Y为每天因为潮汐引起的水位高度
数字化后的信号,每一个点代表采样点

Aliasing

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

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

如果我们对它采样,大概会得到。。。这样:

用一条线把采样点连接起来,采集的到的数字波形可以看出明显的上下振动。由于高频噪声和原来的低频真实信号相叠加,最后采集出来的数据和原来的数据相比很难看。这是因为采样的频率远低于噪声信号的频率,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))

输出结果:

ans =

    0.5649

2. 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
产生2Hz的正弦波

下面准备开始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
fft并取绝对值的plot结果

如果仔细观察,发现定点的值是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 / samples(points) in FFT
频率分辨率 = 采样率 / 一次输入FFT转换的采样点数

在我们的例子中, binwidth(频率分辨率) 就是1Hz.

频率分辨率是啥意思呢? 顾名思义,就是FFT转换后频谱图上的X轴(频率)的最小细分单位。或者说频谱图上X轴的最小间隔。

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)



3. 数字滤波器

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

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

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

FIR filter

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

3点滑动平均滤波器

在这个滤波器中,可以看到每次把前三个数据进行平均(分别乘以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');
原始信号和经过moving average filter滤波后的fft频谱

可以看出,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)...

一个三阶IIR滤波器的框图,Z-1可以理解成暂存上一次的一个结果

4. 自动信号检测

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

能量检测器

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

  1. 对信号进行滤波,将有用信号提取出来
  2. 对信号进行整流
  3. 对整流过的信号进行包络
  4. 设置阈值,检测有用信号

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

% 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函数来找到正弦信号什么时候结束


5. 图像处理

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

编辑于 08-25