درک فیلتر میانگین متحرک بوسیله پایتون و متلب

فیلتر میانگین متحرک (Moving Average – MA)، یک فیلتر با پاسخ ضربه محدود (FIR) ساده پایین‌گذر (Lowpass) است که معمولا برای هموار کردن (Smoothing) آرایه‌ای از نمونه‌های یک سری داده یا سیگنال استفاده میشود. این فیلتر، $L$ نمونه ورودی را در هر زمان دریافت کرده و میانگین آنها را به عنوان یک نقطه در سیگنال خروجی، محاسبه میکند. این فیلتر، یک ساختار فیلتر پایین‌گذر خیلی ساده است که برای دانشمندان و مهندسین برای حذف اجزاء نویزی ناخواسته از داده مورد نظر، بسیار کارآمد است.

با افزایش طول فیلتر (پارامتر $L$ )، میزان صافی و هموار بودن خروجی، افزایش می‌یابد، و از طرفی، تغییرات ناگهانی و شدید در داده نیز واضح‌تر میشوند. از این موضوع چنین میتوان استنباط کرد که این فیلتر پاسخ حوزه زمان خیلی خوبی دارد اما پاسخ فرکانسی آن، ضعیف است.

فیلتر MA، سه کار مهم انجام میدهد:

– $L$ نمونه ورودی را دریافت کرده، میانگین آنها را محاسبه میکند و نتیجه را به صورت یک نقطه در خروجی، تولید میکند.

– به دلیل محاسبات انجام شده در فرآیند فیلتر کردن، مقدار مشخصی از تاخیر در خروجی نسبت به ورودی، ایجاد میشود.

– فیلتر به شکل یک فیلتر پایین‌گذر (با پاسخ حوزه فرکانس ضعیف و پاسخ حوزه زمان خوب)، عمل میکند.

پیاده‌سازی

معادله تفاضلی (Difference) برای فیلتر میانگین متحرک $L$ -نقطه‌ای که ورودی آن با بردار $x$ و خروجی میانگین گرفته شده با $y$ ، نمایش داده میشوند، به شکل زیر است:

$y[n] = \frac{1}{L} \sum_{k=0}^{L-1} x[n-k]$

برای مثال، یک فیلتر FIR میانگین متحرک 5-نقطه‌ای، نمونه فعلی و 4 نمونه قبلی ورودی را گرفته و مقدار میانگین آنها را محاسبه میکند. این عمل که در معادله زیر برای رابطه بین ورودی و خروجی در حوزه زمان گسسته، آمده است، در شکل زیر نیز نمایش داده شده است:

$y[n] = \frac{1}{5}(x[n]+x[n-1]+x[n-2]+x[n-3]+x[n-4])$

شکل 1- فیلتر FIR میانگین متحرک 5-نقطه‌ای

تاخیر واحد در شکل بالا، توسط یکی از دو راه‌حل زیر، قابل تحقق است:

– نمایش نمونه‌های ورودی به شکل یک آرایه در حافظه رایانه و در نهایت پردازش آنها

– استفاده از شیفت رجیسترهای (shift register) با ساختار فلیپ فلاپی D-تایی (D-Flip flop) برای پیاده‌سازی سخت‌افزاری. اگر هر مقدار گسسته ورودی $x$ ، به صورت یک داده 12 بیتی از مبدل آنالوگ به دیجیتال بدست آید، ما نیاز به 4 مجموعه از فلیپ فلاپهای 12-تایی داریم تا بتوانیم فیلتر میانگین متحرک 5-نقطه‌ای در شکل بالا را پیاده‌سازی کنیم.

تبدیل Z و تابع تبدیل

در پردازش سیگنال، تاخیر در یک سیگنال $x[n]$ به میزان $k$ نمونه، معادل ضرب تبدیل Z در عبارت $z^{-k}$ است. با علم به این موضوع، میتوانیم تبدیل Z معادله فیلتر میانگین متحرک 5-نقطه‌ای که در بالا آمده است را به شکل زیر بدست آورد:

$\begin{aligned} y[n] = 0.2(x[n]+x[n-1]+x[n-2]+x[n-3]+x[n-4]) \\ Y[z]=0.2(z^0 + z^{-1}+z^{-2}+z^{-3}+z^{-4})X(z) \\ Y[z]=0.2(1 + z^{-1}+z^{-2}+z^{-3}+z^{-4})X(z) \end{aligned}$

به طور مشابه، تبدیل Z فرم عمومی معادله فیلتر میانگین متحرک $L$ -نقطه‌ای، به شکل زیر حاصل میشود:

$Y(z) = \left(\frac{1}{L} \sum_{k=0}^{L-1} z^{-k} \right) X(z)$

تابع تبدیل، رابطه ورودی-خروجی سیستم را برای فیلتر میانگین متحرک $L$ -نقطه‌ای، مشخص میکند که به شکل زیر بدست می‌آید:

$H(z) = \frac{Y(z)}{X(z)} = \frac{1}{L}\sum_{k=0}^{L-1}z^{-k}$

شبیه‌سازی فیلتر در متلب و پایتون

در متلب (MATLAB)، میتوانیم از تابع filter یا conv (کانولوشن)، برای پیاده‌سازی فیلتر FIR میانگین متحرک استفاده کنیم.

به طور کلی، تبدیل Z فیلتر خروجی فیلتر زمان-گسسته یا همان $Y(z)$ ، با تبدیل Z ورودی $x[n]$ ، یا $X(z)$ به شکل زیر مرتبط میشود:

$H(z) = \frac{Y(z)}{X(z)} = \frac{b_0+b_1 z^{-1}+b_2 z^{-2}+...+b_n z^{-n}}{a_0+a_1 z^{-1}+a_2 z^{-2}+...+a_m z^{-m}}$

که در آن ${b_i}$ و ${a_i}$ ، ضرایب فیلتر هستند و مرتبه فیلتر هم مقدار ماکزیمم $m$ و $n$ است.

برای پیاده‌سازی معادله بالا با استفاده از تابع filter در متلب، به شکل زیر عمل میکنیم:

B = [b_0, b_1, b_2, ..., b_n] %numerator coefficients
A = [a_0, a_1, a_2, ..., a_m] %denominator coefficients
y = filter(B,A,x) %filter input x and get result in y

از معادله تفاضلی و تابع تبدیل فیلتر میانگین متحرک $L$ -نقطه‌ای، میتوان متوجه شد که مقدار ضرایب صورت، ${b_i}$ و ضرایب مخرج، ${a_i}$ برابر است با:

${b_i}={1/L, 1/L, ..., 1/L} \ \ u=0, 1, ..., L-1 \ \ \ \ \ \ \ \ \ \ {a_i} = 1 \ \ i=0$

بنابراین، فیلتر میانگین متحرک 5-نقطه‌ای به شکل زیر میتواند کدنویسی شود:

B = [0.2, 0.2, 0.2, 0.2, 0.2] %numerator coefficients
A = [1] %denominator coefficients
y = filter(B,A,x) %filter input x and get result in y

ضرایب صورت در تابع تبدیل فیلتر میانگین متحرک، به شکل مناسب‌تری به شرح زیر قابل بیان است:

L = 5
B = ones(1,L)/L %numerator coefficients
A = [1] %denominator coefficients
x = rand(1,10) %random samples for x
y = filter(B,A,x) %filter input x and get result in y

اگر از تابع conv برای پیاده‌سازی فیلتر میانگین متحرک استفاده کنیم، کد زیر باید نوشته شود:

L = 5;
x = rand(1,10) %random samples for x;
y = conv(x,ones(1,L)/L)

در پیاده‌سازی به شکل یک فیلتر FIR، بین تابع conv و filter یک اختلاف مهم وجود دارد. تابع conv، حاصل کانولوشن کامل را بدست میدهد که طول آن برابر با length(x)+L-1 است. درحالیکه، تابع filter خروجی با طول برابر با طول سیگنال ورودی را میدهد.

در پایتون، عمل فیلتر کردن، توسط تابع lfilter و convolve که از طریق کتابخانه پردازش سیگنال scipy در دسترس است، انجام میشود. کد معادل پایتون برای فیلتر میانگین متحرک به شرح زیر است:

import numpy as np
from scipy import signal
L=5 #L-point filter
b = (np.ones(L))/L #numerator co-effs of filter transfer function
a = np.ones(1)  #denominator co-effs of filter transfer function
x = np.random.randn(10) #10 random samples for x
y = signal.convolve(x,b) #filter output using convolution

y = signal.lfilter(b,a,x) #filter output using lfilter function

نمودار صفر-قطب و پاسخ فرکانسی

نمودار صفر-قطب یک تابع تبدیل فیلتر $H(z)$ ، مکان صفرها و قطبهای فیلتر را در صفحه z نمایش میدهد. در نمودار صفر-قطب، صفرها در مکانهایی (فرکانسهایی) اتفاق می‌افتند که $H(z) \to 0$ و قطبها در مکانهایی (فرکانسهایی) اتفاق می‌افتند که $H(z) \to \infty$ . به طور معادل، صفرها در فرکانسهایی اتفاق می‌افتند که صورت تابع تبدیل صفر میشود و قطبها در فرکانسهایی اتفاق می‌افتند که مخرج تابع تبدیل صفر میشود.

در نمودار صفر-قطب، مکان قطبها معمولا توسط علامت $\times$ و صفرها توسط علامت $o$ نمایش داده میشوند. قطبها و صفرهای یک تابع تبدیل، به طور موثر، پاسخ سیستم را تعیین کرده و پایداری و عملکرد سیستم فیلترینگ را مشخص میکنند.

در متلب، نمودار صفر-قطب و پاسخ فرکانسی یک فیلتر میانگین متحرک $L$ -نقطه‌ای به شکل زیر قابل دستیابی است:

L=11
zplane([ones(1,L)]/L,1); %pole-zero plot
w = -pi:(pi/100):pi; %to plot frequency response
freqz([ones(1,L)]/L,1,w); %plot magnitude and phase response

شکل 2- نمودار صفر-قطب فیلتر میانگین متحرک 11-نقطه‌ای

شکل 3- دامنه و فاز پاسخ فرکانسی فیلتر میانگین متحرک 11-نقطه‌ای

برای نسخه پایتون کد بالا، داریم:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

L=11 #L-point filter
b = (np.ones(L))/L #numerator co-effs of filter transfer function
a = np.ones(1)  #denominator co-effs of filter transfer function

w, h = signal.freqz(b,a)
plt.subplot(2, 1, 1)
plt.plot(w, 20 * np.log10(abs(h)))
plt.ylabel('Magnitude [dB]')
plt.xlabel('Frequency [rad/sample]')

plt.subplot(2, 1, 2)
angles = np.unwrap(np.angle(h))
plt.plot(w, angles)
plt.ylabel('Angle (radians)')
plt.xlabel('Frequency [rad/sample]')
plt.show()

شکل 4- پاسخ ضربه، نموار صفر-قطب و اندازه و دامنه پاسخ فرکانسی فیلتر میانگین متحرک 11-نقطه‌ای

نمونه کاربرد فیلتر

دو شکل زیر به ترتیب سیگنال ورودی و خروجی حوزه زمان فیلتر میانگین متحرک $L$ -نقطه‌ای و پاسخ فرکانسی آن را برای طولهای مختلف نشان میدهند. در شکل اول، یک سیگنال موج مربعی نویزی از فیلتر عبور داده شده و پاسخ حوزه زمان آن نمایش داده شده است.

شکل 5- پردازش یک سیگنال توسط فیلتر میانگین متحرک با طولهای مختلف

در تصویر بالا، نمودار اول، سیگنال ورودی نویزی است و هدف ما کاهش نویز به میزان حداکثری است. نمودار بعدی، خروجی فیلتر میانگین متحرک 3-نقطه‌ای است. از این شکل میتوان چنین استنباط کرد که این فیلتر 3-نقطه‌ای چندان در حذف نویز موفق نبوده است. تعداد نودهای فیلتر را به 10 افزایش میدهیم و با توجه به شکل بعدی، واضح است که نویز در خروجی به شکل قابل توجهی کاهش یافته است.

با طول $L=51$ ، (در شکل بعدی) نویز تقریبا به صفر میرسد ولی بازه زمانی تغییر سطح سیگنال موج مربعی نیز به شدت افزایش یافته است (شیب خط بین دو سطح سیگنال را با تغییرات ایده‌آل در سیگنال ورودی برای بازه مشابه، مقایسه کنید).

شکل 6- پاسخ فرکانسی فیلتر میانگین متحرک با طولهای مختلف

با توجه به پاسخ فرکانسی فیلتر با طولهای کمتر ( $L=3, L=10$ )، میتوان این موضوع را تایید کرد که تغییر سطح باند عبور به باند توقف بسیار کند انجام میشود و تضعیف باند توقف نیز چندان مناسب نیست. با توجه به این باند توقف، واضح است که فیلتر میانگین متحرک در جدا کردن فرکانسها از یکدیگر نمیتواند موفق عمل کند. با افزایش طول فیلتر به 51 نیز تغییرات حوزه زمان قابل کنترل نیست. بنابراین، عملکرد خوب در حوزه زمان منجر به عملکرد ضعیف در حوزه فرکانس میشود و بالعکس. به همین دلیل همواره با توجه به شرایط و نیاز سیستم، طراحی بهینه انجام میشود.

منبع: www.gaussianwaves.com

درک فیلتر میانگین متحرک بوسیله پایتون و متلب

دیدگاه ها (0)

لغو پاسخ