《电子技术应用》
您所在的位置:首页 > 微波|射频 > 解决方案 > Matlab 数字滤波入门

Matlab 数字滤波入门

2019-07-23

  模拟与数字信号

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

微信图片_20190724123632.jpg

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

微信图片_20190724123634.jpg

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

  Aliasing

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

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

微信图片_20190724123731.jpg

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

微信图片_20190724123734.jpg

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

微信图片_20190724123737.jpg

  避免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轴坐标从频率变成时间

微信图片_20190724123924.jpg

  先用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

微信图片_20190724123958.jpg

  如果仔细观察,发现定点的值是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

微信图片_20190724124054.jpg

  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

微信图片_20190724124057.jpg

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

微信图片_20190724124100.jpg

  好了,大功告成

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

  之前说过,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滤波器的一种。在一些应用中,有一个窗口,每一次新的数据进来都在窗口进行一次平局然后输出。

微信图片_20190724124104.jpg

  在这个滤波器中,可以看到每次把前三个数据进行平均(分别乘以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')

微信图片_20190724124300.jpg

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

微信图片_20190724124304.jpg

  可以看出,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指的是后馈系数。

微信图片_20190724124430.jpg

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

  FDAtool

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

  IIR 滤波器原理

微信图片_20190724124528.jpg

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

微信图片_20190724124531.jpg

  自动信号检测

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

  能量检测器

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

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

  对信号进行整流

  对整流过的信号进行包络

  设置阈值,检测有用信号

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

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

微信图片_20190724124621.jpg

  这个图把所有每一步的波形全部显示出来

  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)

微信图片_20190724124651.jpg

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

  图像处理

微信图片_20190724124722.jpg

  图形实际上可以看做三重矩阵,因为每个像素是RGB三个颜色分量,每个像素都用三个值描述,所以是三重矩阵。让我们先搞一张最简单的,使用向量建立起来的图像:

微信图片_20190724124726.jpg

  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的列组成的

微信图片_20190724124803.jpg

  把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(离散傅里叶变换)

微信图片_20190724124826.jpg

  离散傅里叶变换公式

  其中

  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输出的第一个点为:

微信图片_20190724124828.jpg

  计算FFT的第一个点

微信图片_20190724124831.jpg

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

  参考

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