《電子技術應用》
您所在的位置:首頁 > 微波|射頻 > 解決方案 > 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


重庆时时是合法的吗