فیلترها
در فصلِ حوزهٔ زمان دیدیم که کانولوشن میتواند سیگنال را هموار کند. این، نمونهای ساده از یک فیلتر (filter) بود. بهطورِ کلی، فیلتر ابزاری است که برخی بسامدها را از سیگنال عبور میدهد و برخی دیگر را تضعیف یا حذف میکند. فیلترها در علوم اعصاب نقشِ مرکزی دارند: برای حذفِ نوفهٔ خطِ برق (۵۰ یا ۶۰ هرتز)، برای جداکردنِ باندهای مغزی (آلفا، بتا، گاما)، و برای حذفِ روندِ آهستهٔ پسزمینه از ثبتها.
انواع فیلتر
بر پایهٔ اینکه کدام بسامدها را عبور میدهند، فیلترها را به چند دسته تقسیم میکنیم:
- فیلترِ پایینگذر (low-pass): بسامدهای پایینتر از یک بسامدِ مرزی (cutoff) را عبور میدهد و بسامدهای بالا را تضعیف میکند. برای هموارسازی و حذفِ نوفهٔ پربسامد به کار میرود.
- فیلترِ بالاگذر (high-pass): برعکس، بسامدهای بالا را عبور میدهد و بسامدهای پایین (مثلاً روندِ آهسته) را حذف میکند.
- فیلترِ میانگذر (band-pass): تنها بسامدهای میانِ دو مرز را عبور میدهد. برای جداکردنِ یک باندِ بسامدیِ خاص (مثلاً باندِ آلفا) آرمانی است.
- فیلترِ میاننگذر (band-stop یا notch): برعکسِ میانگذر، تنها یک باندِ باریک را حذف میکند. کاربردِ کلاسیکِ آن، حذفِ نوفهٔ ۵۰ هرتزیِ خطِ برق است.
پاسخ بسامدی
رفتارِ یک فیلتر را با پاسخِ بسامدیِ آن توصیف میکنیم: نموداری که نشان میدهد فیلتر به هر بسامد چه بهرهای (gain) میدهد. بهرهٔ نزدیک به ۱ یعنی آن بسامد تقریباً دستنخورده عبور میکند، و بهرهٔ نزدیک به ۰ یعنی آن بسامد حذف میشود. ناحیهای که فیلتر عبور میدهد باندِ عبور (passband) و ناحیهای که حذف میکند باندِ توقف (stopband) نام دارد.
در عمل، گذار از باندِ عبور به باندِ توقف هرگز کاملاً تند نیست؛ همیشه یک ناحیهٔ گذارِ تدریجی وجود دارد. مرتبهٔ فیلتر تعیین میکند که این گذار چقدر تند است: مرتبهٔ بالاتر، گذارِ تندتر، اما به بهای پیچیدگیِ بیشتر و احتمالِ ناپایداری.
فیلترهای FIR و IIR
پیش از آنکه به ساختِ فیلترها بپردازیم، باید بدانیم هر فیلترِ خطی به یکی از دو خانوادهٔ بنیادی تعلق دارد. کلیدِ تمایز، پاسخِ ضربهای (impulse response) است: خروجیِ فیلتر وقتی ورودی یک ضربهٔ تنها \((1, 0, 0, \dots)\) باشد. پاسخِ ضربهای، اثرِ انگشتِ فیلتر است؛ هر فیلترِ خطی با کانوالوکردنِ ورودی با پاسخِ ضربهایِ خود کار میکند.
پاسخِ ضربهای چیست و چرا همهچیز را تعیین میکند؟
ضربه (impulse) سادهترین سیگنالِ ممکن است: تنها در لحظهٔ صفر مقدارِ ۱ دارد و در بقیهٔ جاها صفر است، یعنی \(\delta[n] = (1, 0, 0, 0, \dots)\). در زمانِ گسسته به آن واحدِ ضربه یا دلتای کرونکر (Kronecker delta) میگویند؛ همتای پیوستهاش دلتای دیراک (Dirac delta) \(\delta(t)\) است—که بهجای یک دنباله، یک تابعِ تعمیمیافته با پهنای صفر و سطحِ زیرِ منحنیِ واحد است. پاسخِ ضربهای \(h[n]\)، خروجیِ فیلتر در پاسخ به همین ورودیِ ضربه است.
چرا این یک سیگنالِ کوچک، کلِ رفتارِ فیلتر را تعیین میکند؟ نکته اینجاست که هر سیگنال را میتوان بهصورتِ مجموعی از ضربههای جابهجاشده و مقیاسخورده نوشت:
حال اگر فیلتر خطی و ناوردا به انتقال (LTI) باشد، دو خاصیت داریم: پاسخ به مجموع، مجموعِ پاسخهاست (خطیبودن)، و پاسخ به ضربهٔ جابهجاشده، همان پاسخِ ضربهایِ جابهجاشده است (ناوردایی به انتقال). با ترکیبِ این دو، خروجی برای هر ورودی چنین میشود:
یعنی دانستنِ \(h\) بهتنهایی، خروجی را برای هر ورودیِ ممکن میدهد—و آن عمل، دقیقاً همان کانولوشن است. به همین دلیل \(h\) را «اثرِ انگشتِ» فیلتر مینامیم: همان هستهٔ کانولوشن است، و تبدیلِ فوریهاش همان پاسخِ بسامدیِ فیلتر.
تفاوتِ FIR و IIR هم در همینجا ریشه دارد: در فیلترِ FIR پاسخِ ضربهای پس از چند نمونه دقیقاً صفر میشود (متناهی)، اما در فیلترِ IIR بهخاطرِ بازخورد، تا بینهایت ادامه مییابد و تنها بهتدریج میرا میشود (نامتناهی):
import numpy as np
from scipy import signal as sig
impulse = np.zeros(60); impulse[0] = 1.0 # the impulse (1, 0, 0, ...)
h_fir = sig.lfilter(np.ones(11)/11, [1], impulse) # FIR: finite
b, a = sig.butter(4, 0.12)
h_iir = sig.lfilter(b, a, impulse) # IIR: infinite (decaying)
پاسخِ ضربهای و تابعِ گرین (Green's function)
آنچه در پردازشِ سیگنال پاسخِ ضربهای مینامیم، در فیزیک و حلِّ معادلاتِ دیفرانسیل تابعِ گرین (Green's function) نام دارد—همان شیء، با دو نام. تابعِ گرینِ یک عملگرِ خطیِ \(L\)، پاسخِ آن به یک چشمهٔ دلتای دیراک است:
و پاسخ به هر ورودیِ دلخواه از اصلِ برهمنهی بهدست میآید: \(u(x) = \int G(x, y)\, f(y)\, dy\). برای سامانهٔ LTI این انتگرال دقیقاً همان کانولوشنِ \(u = f * h\) میشود. دلیلش هم همان است که پاسخِ ضربهای کلِ فیلتر را تعیین میکرد: هر ورودی مجموعی از ضربههای جابهجاشده است، و خروجی مجموعِ پاسخهای ضربهایِ (تابعهای گرینِ) جابهجاشده—که تعریفِ کانولوشن است:
این پیوند چنان ژرف است که سزاوارِ نگاهی مفصلتر است: در بخشِ پایانیِ فصل، «تابعِ گرین: نگاهی ژرفتر»، آن را در معادلهٔ گرما، الکترواستاتیک، شرایطِ مرزی، کدِ گسسته، حوزهٔ بسامد، و حتی هوشِ مصنوعی دنبال میکنیم.
فیلترِ FIR (پاسخِ ضربهایِ متناهی، Finite Impulse Response): پاسخِ ضربهای طولِ متناهی دارد. خروجی تنها ترکیبی خطی از مقادیرِ گذشته و حالِ ورودی است:
این دقیقاً همان کانولوشنی است که در فصلِ حوزهٔ زمان دیدیم: ضرایبِ \(b_k\) همان هستهٔ کانولوشن (پاسخِ ضربهای) هستند. میانگینِ متحرکی که در آن فصل ساختیم، سادهترین فیلترِ FIR است (هستهای جعبهای).
فیلترِ IIR (پاسخِ ضربهایِ نامتناهی، Infinite Impulse Response): پاسخِ ضربهای طولِ نامتناهی دارد و با یک معادلهٔ تفاضلی توصیف میشود که علاوه بر ورودی، به مقادیرِ گذشتهٔ خروجی نیز وابسته است (جملهٔ بازخورد):
همین بازخورد است که IIR را نیرومند (با ضرایبِ کم، گذارِ تند) اما در عوض پیچیدهتر و مستعدِ ناپایداری میکند. وقتی در ادامه scipy.signal.butter را فرامیخوانیم، همین دو دستهٔ ضرایب—\(b\) (پیشخور) و \(a\) (بازخورد)—را برمیگرداند، و filtfilt آنها را با همین معادلهٔ تفاضلی اعمال میکند.
برای دیدنِ تفاوت، یک فیلترِ پایینگذرِ FIR (با تابعِ firwin) و یک فیلترِ پایینگذرِ IIR باترورث را با مرزِ یکسانِ ۳۰ هرتز میسازیم و مقایسه میکنیم:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as sig
fs = 500.0
fc = 30.0
b_fir = sig.firwin(61, fc, fs=fs) # FIR: 61 taps = the impulse response
b_iir, a_iir = sig.butter(4, fc, fs=fs) # IIR: Butterworth, order 4
w1, h1 = sig.freqz(b_fir, [1], fs=fs, worN=2000)
w2, h2 = sig.freqz(b_iir, a_iir, fs=fs, worN=2000)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3.6))
ax1.stem(np.arange(len(b_fir)), b_fir, linefmt="tab:blue", markerfmt=" ", basefmt=" ")
ax1.set_title("FIR impulse response (61 taps)")
ax1.set_xlabel("tap k"); ax1.set_ylabel("b[k]")
ax2.plot(w1, np.abs(h1), color="tab:blue", label="FIR (61 taps)")
ax2.plot(w2, np.abs(h2), color="tab:red", label="IIR Butterworth (order 4)")
ax2.axvline(fc, color="gray", ls=":", lw=1)
ax2.set_xlim(0, 100); ax2.set_xlabel("frequency (Hz)"); ax2.set_ylabel("gain")
ax2.set_title("frequency response"); ax2.legend()
plt.tight_layout()
plt.show()
کدام را برگزینیم؟
فیلترِ FIR همیشه پایدار است و میتواند فازِ خطی داشته باشد (تأخیرِ یکسان برای همهٔ بسامدها)، اما برای گذارِ تند به ضرایبِ بسیار زیاد نیاز دارد. فیلترِ IIR با ضرایبِ بسیار کمتر به گذارِ تند میرسد، اما فازش خطی نیست و میتواند ناپایدار شود. در علوم اعصاب، وقتی زمانبندیِ دقیق مهم است، فازِ خطیِ FIR ارزشمند است؛ وقتی کارایی مهم است، IIR (مانندِ باترورث) رایجتر است.
کانولوشن و فیلتر: دو روی یک سکه
در فصلِ حوزهٔ زمان، کانولوشن و میانگینِ متحرک را ساختیم؛ در این فصل از فیلترهای پایینگذر، بالاگذر و میانگذر سخن میگوییم. آیا اینها دو چیزِ متفاوتاند؟ پاسخ، که در بخشِ پیش هم به آن اشاره شد، این است: نه—اینها دو نگاه به یک عملِ واحدند.
شباهت. هر فیلترِ خطی، در دلِ خود یک کانولوشن با پاسخِ ضربهای (هسته) است، و برعکس، هر کانولوشن با یک هستهٔ ثابت، یک فیلترِ خطی است. پس «میانگینِ متحرک» و «فیلترِ پایینگذر» میتوانند یک چیز باشند.
تفاوت—در واژگان، نه در ماهیت. آنچه فرق میکند، زاویهٔ نگاه است:
- وقتی فیلتر را با هستهاش نام میبریم—«فیلترِ میانگینِ متحرک»، «فیلترِ گاوسی»—داریم آن را از منظرِ حوزهٔ زمان توصیف میکنیم: با چه چیزی کانوالو میکنیم.
- وقتی آن را «پایینگذر»، «بالاگذر» یا «میانگذر» مینامیم، داریم آن را از منظرِ حوزهٔ بسامد توصیف میکنیم: کدام بسامدها را عبور میدهد.
پلِ میانِ این دو نگاه، قضیهٔ کانولوشن است: تبدیلِ فوریهٔ هسته، همان پاسخِ بسامدیِ فیلتر است. پس شکلِ هسته در حوزهٔ زمان، نوعِ فیلتر را در حوزهٔ بسامد تعیین میکند. بیایید چند هسته را کنارِ پاسخِ بسامدیِ آنها ببینیم:
import numpy as np
import matplotlib.pyplot as plt
fs = 1000.0
def freq_response(kernel, nfft=4096):
H = np.fft.rfft(kernel, nfft)
f = np.fft.rfftfreq(nfft, 1/fs)
mag = np.abs(H)
return f, mag/mag.max()
box = np.ones(11)/11 # moving average
xg = np.linspace(-3, 3, 21); g = np.exp(-xg**2/2); gauss = g/g.sum() # Gaussian
diff = np.array([1.0, -1.0]) # first difference
xx = np.linspace(-4, 4, 41)
gn = np.exp(-xx**2/(2*0.5**2)); gn /= gn.sum()
gw = np.exp(-xx**2/(2*1.5**2)); gw /= gw.sum()
dog = gn - gw # difference of Gaussians
rows = [("moving average (box)", "low-pass", box),
("Gaussian", "low-pass (smooth)", gauss),
("first difference", "high-pass", diff),
("difference of Gaussians", "band-pass", dog)]
fig, axes = plt.subplots(4, 2, figsize=(10, 8))
for r, (name, ftype, k) in enumerate(rows):
axes[r, 0].stem(np.arange(len(k))-len(k)//2, k, linefmt="tab:blue",
markerfmt=" ", basefmt=" ")
axes[r, 0].set_title(f"kernel: {name}", fontsize=10); axes[r, 0].set_ylabel("h[k]")
f, H = freq_response(k)
axes[r, 1].plot(f, H, color="tab:red"); axes[r, 1].set_xlim(0, fs/2)
axes[r, 1].set_title(f"frequency response -> {ftype}", fontsize=10)
axes[r, 1].set_ylabel("gain")
plt.tight_layout()
plt.show()
این تصویر همهٔ مفاهیم را به هم گره میزند: هستهٔ جعبهای که در فصلِ پیش برای هموارسازی به کار بردیم، یک فیلترِ پایینگذر بوده است—و موجهای جانبیِ آن (همان سینوسیشکلِ ناشی از لبههای تیزِ جعبه) دلیلِ آن است که هستهٔ گاوسی برای هموارسازی بهتر است. هستهٔ تفاضلی (که تغییراتِ تند را برجسته میکند) یک فیلترِ بالاگذر است، و تفاضلِ دو گاوسی یک فیلترِ میانگذر.
کجا این همارزی میشکند؟
گزارهٔ «فیلتر = کانولوشن با یک هستهٔ متناهی» دقیقاً برای فیلترهای FIR برقرار است. برای فیلترهای IIR (مانندِ باترورث در بخشِ بعد)، پاسخِ ضربهای نامتناهی است؛ نمیتوان با یک هستهٔ متناهی کانوالو کرد، و بهجای آن از معادلهٔ تفاضلیِ بازگشتی استفاده میکنیم. پس دقیقترین بیان این است: هر فیلترِ خطی کانولوشن با پاسخِ ضربهایِ خود است—که برای FIR متناهی و مستقیماً قابلاستفاده، و برای IIR نامتناهی است.
پیوند با مدلهای سریهای زمانی: ARIMA و SARIMA
همین آجرهای سازنده، در آمارِ سریهای زمانی نیز—با نامهایی دیگر—دوباره ظاهر میشوند. مدلهای ARIMA برای مدلسازی و پیشبینیِ یک سری به کار میروند (نه برای فیلتر کردنِ آن)، اما ماشینِ ریاضیِ پشتشان همان است:
- AR (خودبازگشتی، autoregressive): \(y_n = \sum_{i=1}^{p} \varphi_i\, y_{n-i} + \varepsilon_n\). این دقیقاً همان جملهٔ بازخوردیِ یک فیلترِ IIR است (وابستگی به مقادیرِ گذشتهٔ خروجی)، اما اینبار با ورودیِ نوفهٔ سفید \(\varepsilon_n\). به بیانِ دیگر، یک فرایندِ AR، خروجیِ یک فیلترِ IIR است که با نوفه رانده میشود.
- MA (میانگینِ متحرک—اما نه آن میانگینِ متحرک!): در ARIMA، \(y_n = \varepsilon_n + \sum_{i=1}^{q} \theta_i\, \varepsilon_{n-i}\)، یعنی ترکیبی وزندار از خطاهای گذشته، نه میانگینی از خودِ سیگنال. این یک ساختارِ پیشخور (FIR) روی نوفه است و با فیلترِ هموارسازِ «میانگینِ متحرک»ی که در این فصل ساختیم متفاوت است؛ تنها همناماند.
- ARMA = ترکیبِ هر دو = همان معادلهٔ تفاضلیِ کاملِ یک فیلتر (هم قطبها، هم صفرها)، که با نوفه رانده میشود.
- I (انباشته، integrated): ARIMA پیش از مدلسازی، سری را \(d\) بار تفاضل میگیرد تا روندِ ناایستا حذف شود. تفاضلِ نخست (\(y_n - y_{n-1}\)) دقیقاً همان هستهٔ بالاگذرِ \([1, -1]\) است که در گالریِ بالا دیدیم: حذفِ روند، نوعی فیلتر کردنِ بالاگذر است.
- SARIMA: نسخهٔ فصلیِ ARIMA که جملهها و تفاضلگیری را در یک دورهٔ فصلیِ \(s\) نیز میافزاید. تفاضلِ فصلی (\(y_n - y_{n-s}\)) یک فیلترِ شانهای (comb) است که چرخهٔ فصلی و همنواهایش را حذف میکند.
تفاوتِ اصلی در هدف است، نه در ابزار: فیلترها معمولاً برای پاکسازی و استخراج به کار میروند، و ARIMA/SARIMA برای مدلسازی و پیشبینی. تعیینِ مرتبههای این مدلها (با کمکِ خودهمبستگی و خودهمبستگیِ جزئی) و مفهومِ ایستایی، موضوعِ فصلِ تحلیلِ سریهای زمانی است.
فیلترهای باترورث در scipy
یکی از پرکاربردترین خانوادههای فیلتر، باترورث (Butterworth) است که در باندِ عبور پاسخی بسیار هموار (بدونِ موج) دارد. کتابخانهٔ scipy.signal ساختنِ این فیلترها را ساده میکند. تابعِ butter ضرایبِ فیلتر را میسازد و filtfilt آن را بر سیگنال اعمال میکند (دوبار، یکبار رو به جلو و یکبار رو به عقب، تا فاز جابهجا نشود).
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as sig
def butter_filter(x, fs, cutoff, btype, order=4):
# design a Butterworth filter and apply it with zero phase shift
b, a = sig.butter(order, cutoff, btype=btype, fs=fs)
return sig.filtfilt(b, a, x)
# a signal with three components: 2 Hz, 15 Hz and 50 Hz
fs = 500.0
t = np.arange(0, 2, 1/fs)
x = (np.sin(2*np.pi*2*t)
+ 0.7*np.sin(2*np.pi*15*t)
+ 0.5*np.sin(2*np.pi*50*t))
low = butter_filter(x, fs, 5, "low") # keeps the 2 Hz component
band = butter_filter(x, fs, [8, 25], "band") # keeps the 15 Hz component
high = butter_filter(x, fs, 30, "high") # keeps the 50 Hz component
fig, axes = plt.subplots(4, 1, figsize=(8.5, 7), sharex=True)
axes[0].plot(t, x, color="gray"); axes[0].set_title("original: 2 + 15 + 50 Hz")
axes[1].plot(t, low, color="tab:blue"); axes[1].set_title("low-pass: keeps 2 Hz")
axes[2].plot(t, band, color="tab:green"); axes[2].set_title("band-pass: keeps 15 Hz")
axes[3].plot(t, high, color="tab:red"); axes[3].set_title("high-pass: keeps 50 Hz")
axes[3].set_xlabel("time t (s)"); axes[3].set_xlim(0, 1)
plt.tight_layout()
plt.show()
میتوان پاسخِ بسامدیِ این فیلترها را نیز مستقیماً رسم کرد تا ببینیم هر کدام چه بسامدهایی را عبور میدهند. تابعِ freqz پاسخِ بسامدی را میدهد:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as sig
fs = 500.0
designs = [(5, "low", "low-pass (5 Hz)"),
([8, 25], "band", "band-pass (8-25 Hz)"),
(30, "high", "high-pass (30 Hz)")]
for cutoff, btype, label in designs:
b, a = sig.butter(4, cutoff, btype=btype, fs=fs)
w, h = sig.freqz(b, a, fs=fs, worN=2000)
plt.plot(w, np.abs(h), label=label)
plt.xlim(0, 80)
plt.xlabel("frequency (Hz)")
plt.ylabel("gain")
plt.legend()
plt.show()
فیلتر کردن، یک کانولوشن است
هر فیلترِ خطی، در دلِ خود یک کانولوشن است (فصلِ حوزهٔ زمان): سیگنال را با پاسخِ ضربهایِ فیلتر کانوالو میکنیم. بهطورِ همارز، در حوزهٔ بسامد، فیلتر کردن یعنی ضربِ طیفِ سیگنال در پاسخِ بسامدیِ فیلتر. این، نمونهٔ دیگری از قضیهٔ کانولوشن است: کانولوشن در حوزهٔ زمان، با ضرب در حوزهٔ بسامد همارز است.
هشدار: فاز و پیچش
اعمالِ فیلتر میتواند فازِ سیگنال را جابهجا کند، یعنی رویدادها را در زمان حرکت دهد. برای تحلیلهایی که زمانبندیِ دقیق مهم است (مانندِ پتانسیلهای وابسته به رویداد در EEG)، از فیلترِ بدونِفاز مانندِ filtfilt استفاده میکنیم که سیگنال را دوبار (جلو و عقب) فیلتر میکند تا جابهجاییِ فاز خنثی شود. همچنین، فیلترِ با مرتبهٔ بسیار بالا ممکن است ناپایدار شود یا در لبههای سیگنال پیچش (artifact) ایجاد کند.
کاربردهای عملی
تا اینجا فیلترها را روی سیگنالهای مصنوعی آزمودیم. حال سه کاربردِ واقعی در پردازشِ دادههای مغزی را میبینیم که هر روز در آزمایشگاه به کار میروند. در هر سه از همان تابعِ butter_filter که بالاتر تعریف کردیم استفاده میکنیم.
فیلترِ پایینگذر برای حذفِ نوفه
پرکاربردترین کاربردِ فیلترِ پایینگذر، حذفِ نوفه است. فرض کنید یک ریتمِ آهستهٔ مغزی (مثلاً ریتمِ تتای ۶ هرتزی) را ثبت کردهایم، اما ثبت با نوفهٔ پربسامدِ فراوانی آلوده است. چون سیگنالِ موردِ علاقهٔ ما آهسته و نوفه تند است، یک فیلترِ پایینگذر میتواند نوفه را بزداید و ریتم را آشکار کند:
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(0)
fs = 1000.0
t = np.arange(0, 2, 1/fs)
clean = np.sin(2*np.pi*6*t) # a 6 Hz rhythm (theta-like)
noisy = clean + 0.8*rng.standard_normal(len(t)) # heavy broadband noise
denoised = butter_filter(noisy, fs, 15, "low") # low-pass at 15 Hz
plt.plot(t, noisy, color="gray", alpha=0.6, lw=0.7, label="noisy")
plt.plot(t, denoised, color="tab:blue", lw=1.6, label="low-pass (15 Hz)")
plt.plot(t, clean, "--", color="tab:red", lw=1.2, label="true 6 Hz")
plt.xlim(0, 1); plt.xlabel("time t (s)"); plt.ylabel("x(t)")
plt.legend()
plt.show()
نکتهٔ کلیدی، انتخابِ بسامدِ مرزی است: باید بهقدری بالا باشد که سیگنالِ موردِ نظر (۶ هرتز) را دستنخورده عبور دهد، و بهقدری پایین که نوفهٔ پربسامد را حذف کند. مرزِ ۱۵ هرتز این بدهبستان را خوب برآورده میکند.
فیلترِ بالاگذر برای یافتنِ اسپایکها
یکی از زیباترین کاربردهای فیلترِ بالاگذر، جداکردنِ اسپایکها (پتانسیلهای عمل) در ثبتهای خارجسلولی است. یک الکترودِ خارجسلولی همزمان دو چیز را ثبت میکند: پتانسیلِ میدانیِ محلی (LFP) که آهسته و بزرگدامنه است، و اسپایکها که تند (پهنای حدودِ ۱ میلیثانیه) و کوچکدامنهاند. در سیگنالِ خام، خطِ پایه با نوسانِ آهستهٔ LFP بالا و پایین میرود، و به همین دلیل نمیتوان با یک آستانهٔ ثابت اسپایکها را یافت. اگر سیگنال را از یک فیلترِ بالاگذر (معمولاً با مرزِ ۳۰۰ هرتز) بگذرانیم، LFP حذف میشود، خطِ پایه صاف میشود، و اسپایکها بهروشنی بیرون میزنند؛ آنگاه یک آستانهٔ ساده آنها را آشکار میکند:
import numpy as np
import matplotlib.pyplot as plt
fs = 30000.0 # 30 kHz, typical extracellular rate
t = np.arange(0, 1.0, 1/fs)
rng = np.random.default_rng(3)
# slow LFP (large, low-frequency) + brief spikes + noise
lfp = 1.8*np.sin(2*np.pi*4*t) + 1.0*np.sin(2*np.pi*9*t)
w = int(0.0015*fs); tt = np.linspace(-1, 1, w)
wav = -np.exp(-(tt*2.2)**2)*tt*3.0; wav = wav/np.max(np.abs(wav)) # ~1.5 ms spike shape
spikes = np.zeros_like(t)
spk_idx = np.sort(rng.choice(np.arange(w, len(t)-w), size=18, replace=False))
for s in spk_idx:
spikes[s:s+w] += 0.5*wav
raw = lfp + spikes + 0.04*rng.standard_normal(len(t))
hp = butter_filter(raw, fs, 300, "high") # high-pass isolates the spikes
# robust threshold (median absolute deviation) and refractory detection
sigma = np.median(np.abs(hp)) / 0.6745
thr = -4*sigma
refractory = int(0.002*fs) # 2 ms dead time
below = hp < thr
candidates = np.where((~below[:-1]) & (below[1:]))[0]
detected = []; last = -refractory
for c in candidates:
if c - last > refractory:
w1 = min(c + int(0.001*fs), len(hp))
detected.append(c + int(np.argmin(hp[c:w1]))) # locate the spike peak
last = c
detected = np.array(detected)
fig, (a1, a2) = plt.subplots(2, 1, figsize=(9, 4.6), sharex=True)
a1.plot(t*1000, raw, color="gray", lw=0.5); a1.set_ylabel("raw")
a1.set_title("raw: a fixed threshold fails — the LFP baseline wanders")
a2.plot(t*1000, hp, color="tab:blue", lw=0.5)
a2.axhline(thr, color="tab:red", ls=":", lw=1, label="threshold")
a2.plot(detected/fs*1000, hp[detected], "v", color="tab:red", ms=6, label="detected")
a2.set_ylabel("high-pass"); a2.set_xlabel("time (ms)")
a2.set_title("high-pass (300 Hz): flat baseline, spikes detected")
a2.legend(); a2.set_xlim(0, 1000)
plt.tight_layout()
plt.show()
آستانهٔ مقاوم
برای تعیینِ آستانه، بهجای انحرافِ معیار از یک برآوردِ مقاوم بهره میبریم: \(\sigma \approx \text{median}(|x|)/0.6745\). این برآورد، برخلافِ انحرافِ معیار، تحتِ تأثیرِ خودِ اسپایکها (که مقادیرِ پرتاند) قرار نمیگیرد و در عملِ مرتبسازیِ اسپایک (spike sorting) استانداردی رایج است. «زمانِ مرده» (refractory) نیز تضمین میکند هر اسپایک تنها یکبار شمرده شود.
فیلترِ میاننگذر برای حذفِ نوفهٔ خطِ برق
نوفهٔ ۵۰ هرتزیِ خطِ برق (۶۰ هرتز در برخی کشورها) آفتِ همیشگیِ ثبتهای الکتروفیزیولوژیک است. چون این نوفه در یک بسامدِ باریکِ مشخص است، یک فیلترِ میاننگذر (notch) آرمانیترین ابزار برای حذفِ آن است: تنها همان باندِ باریک را حذف میکند و بقیهٔ سیگنال را دستنخورده میگذارد. تابعِ scipy.signal.iirnotch این فیلتر را میسازد:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as sig
fs = 1000.0
t = np.arange(0, 3, 1/fs)
brain = np.sin(2*np.pi*10*t) + 0.6*np.sin(2*np.pi*22*t) # 10 + 22 Hz rhythms
line = 1.2*np.sin(2*np.pi*50*t) # 50 Hz mains noise
contaminated = brain + line
b, a = sig.iirnotch(50, Q=30, fs=fs) # notch at 50 Hz
cleaned = sig.filtfilt(b, a, contaminated)
f1, P1 = sig.welch(contaminated, fs=fs, nperseg=1024)
f2, P2 = sig.welch(cleaned, fs=fs, nperseg=1024)
plt.semilogy(f1, P1, color="gray", label="contaminated")
plt.semilogy(f2, P2, color="tab:green", label="notch @ 50 Hz")
plt.xlim(0, 80); plt.xlabel("frequency (Hz)"); plt.ylabel("power")
plt.legend()
plt.show()
تابعِ گرین: نگاهی ژرفتر
در بخشِ فیلترهای FIR و IIR دیدیم که پاسخِ ضربهای و تابعِ گرین یک شیءاند. این پیوند چنان بنیادی است که در سراسرِ فیزیک، مهندسی و یادگیریِ ماشین تکرار میشود. در این بخش آن را با جزئیاتِ بیشتری دنبال میکنیم؛ این مطالب برای فهمِ بقیهٔ کتاب ضروری نیست، اما تصویری یکپارچه از «فیلتر، کانولوشن و معادلهٔ دیفرانسیل» به دست میدهد.
برای دیدنِ پیوند، سه نقش را از هم جدا کنیم: هسته/کرنل (چگونه)—در پردازشِ سیگنال و تصویر، کانولوشن با یک هستهٔ کوچک مقادیرِ همسایه را با هم میآمیزد (مثلاً لغزاندنِ یک ماتریسِ \(3\times3\) روی پیکسلها برای تار کردن یا آشکارسازیِ لبه)؛ تابعِ گرین (چرا)—در فیزیک، تابعِ گرین همان هستهٔ پاسخِ ضربهایِ یک معادلهٔ دیفرانسیلِ خطی است و توصیف میکند که سامانه به یک «ضربهٔ» کامل (یک تلنگرِ نقطهای، یک جهشِ کوتاهِ ولتاژ، یا یک بارِ نقطهای) چگونه پاسخ میدهد؛ و کانولوشن (پیوند)—چون سامانه خطی است، هر ورودیِ پیچیده مجموعی از این ضربههای ریز است، و با کانولوشنِ تابعِ گرین روی ورودی، خروجیِ نهایی (فیلترشده) به دست میآید.
به زبانِ ریاضی، اگر بخواهیم معادلهٔ دیفرانسیلِ خطیِ \(L\,u(x) = f(x)\) را حل کنیم (که \(L\) عملگرِ دیفرانسیلی، \(f\) ورودی، و \(u\) خروجی است)، تابعِ گرین با \(L\,G(x,y) = \delta(x-y)\) تعریف میشود و پاسخِ کامل از انتگرالِ کانولوشن \(u(x) = \int G(x,y)\,f(y)\,dy\) به دست میآید. هم تابعِ گرین و هم فیلترِ خطی، کاملاً به خطی و ناوردا به انتقال (LTI) بودنِ سامانه تکیه دارند: ویژگیهای سامانه با زمان تغییر نمیکنند، و پاسخ به مجموعِ ورودیها، مجموعِ پاسخهاست.
معادلهٔ گرما و تاریِ گاوسی
کلاسیکترین مثال، پخشِ گرما در یک میلهٔ بلند است. اگر در لحظهٔ \(t=0\) توزیعِ دمای دلخواهِ \(f(x)\) را برقرار کنیم، دما از معادلهٔ گرما پیروی میکند:
برای حل با تابعِ گرین، نخست پاسخ به یک چشمهٔ گرمای نقطهای (یک دلتای دیراک در مبدأ) را مییابیم. این پاسخِ بنیادی، تابعِ گرین، یک گاوسی است که با زمان پهنتر میشود:
چون معادله خطی است، پاسخ به هر ورودیِ دلخواه از کانولوشنِ \(f\) با همین تابعِ گرین به دست میآید: \(u(x,t) = (f * G)(x)\). حال اگر در یک زمانِ ثابتِ \(t = T\) عکس بگیریم، \(G\) یک گاوسیِ ثابت در مکان میشود (با \(\sigma^2 = 2kT\))—که دقیقاً فرمولِ فیلترِ گاوسی در پردازشِ سیگنال و تصویر است. پس پخشِ گرما برای مدتِ \(T\)، عیناً همان تاریکردنِ گاوسیِ (هموارسازیِ) ورودیِ اولیه است: نرمافزار وقتی روی تصویری «تاریِ گاوسی» اعمال میکند، در واقع حساب میکند که گرما در یک لحظهٔ کوتاه چگونه میانِ پیکسلها پخش میشود.
این تناظر را میتوان در یک جدول خلاصه کرد:
| مفهوم | فیزیک / معادلهٔ دیفرانسیل | پردازشِ سیگنال و تصویر |
|---|---|---|
| ورودی | پروفایلِ دمای اولیه \(f(x)\) | سیگنال یا تصویرِ نوفهایِ اصلی \(f(x)\) |
| عملگرِ سامانه | معادلهٔ گرما \(\frac{\partial}{\partial t} - k\frac{\partial^2}{\partial x^2}\) | الگوریتمِ تاری/فیلترِ پایینگذر |
| هسته | تابعِ گرین \(G(x,t)\) | پاسخِ ضربهای / هستهٔ گاوسی \(h(x)\) |
| خروجی | دمای پخششده \(u(x,t)\) | سیگنالِ هموار/فیلترشده \(g(x)\) |
الکترواستاتیک: پتانسیلِ کولن
تابعِ گرین تنها مخصوصِ مسائلِ وابسته به زمان نیست. در الکترواستاتیک، سامانهای ایستا داریم که از معادلهٔ پواسون پیروی میکند:
اینجا مشتقِ زمانی نداریم، چون بارها ساکناند. تابعِ گرین \(G(\mathbf{r}, \mathbf{r}')\) پتانسیلِ الکتریکی در نقطهٔ \(\mathbf{r}\) است که از یک بارِ نقطهایِ واحد (یک دلتای دیراک) در \(\mathbf{r}'\) ناشی میشود—یعنی همان پتانسیلِ کولنِ آشنا:
چون نیروهای الکترواستاتیک از برهمنهیِ خطی پیروی میکنند، پتانسیلِ کلِ هر توزیعِ بارِ دلخواهِ \(\rho(\mathbf{r}')\) از کانولوشنِ آن با تابعِ گرین به دست میآید:
این دقیقاً همان مفهومِ معادلهٔ گرماست؛ تنها بهجای ردگیریِ پخشِ گرما در زمان، در حالِ محاسبهٔ گسترشِ ولتاژ در مکان از یک چشمهٔ بار هستیم.
شرایطِ مرزی و روشِ تصویرها
سامانههای واقعی بینهایتبزرگ نیستند؛ دیوار، لبه و مرز دارند، و تابعِ گرین باید شرایطِ مرزی را برآورده کند. مثالِ کلاسیک، یک بارِ نقطهایِ مثبتِ \(+q\) در کنارِ یک صفحهٔ فلزیِ زمینشدهٔ بینهایت است، که در آن \(V = 0\) یک قیدِ مرزیِ سخت است. تابعِ گرینِ فضای آزاد بهتنهایی کار نمیکند، چون دیوارِ فلزی را به حساب نمیآورد. ترفندِ روشِ تصویرها (method of images) چنین است: دیوار را برمیداریم و یک بارِ «آینهای» منفیِ \(-q\) را دقیقاً در نقطهٔ متقارن، آنسوی محلِ دیوار، میگذاریم. تابعِ گرینِ جدید و آگاه از مرز، مجموعِ این دو بار است:
بهطراحی، اگر هر مختصاتی روی صفحهٔ مرزی را در این تابع بگذاریم، دو جمله دقیقاً یکدیگر را خنثی میکنند و \(V = 0\) تضمین میشود.
پیوند با حالتِ مرزی در کد
وقتی این مرزها را به کد یا ماتریس ترجمه میکنیم، باید به برنامه بگوییم لبهها را چگونه مدیریت کند؛ به همین دلیل کتابخانههای فیلتر از ما حالتِ مرزی (boundary mode) میخواهند. دو حالتِ رایج، تناظرِ مستقیم با فیزیک دارند: صفرگذاری (mode='constant', cval=0) فرض میکند بیرونِ مرز همهچیز صفر است—معادلِ همان دیوارِ زمینشده (\(V=0\))؛ و بازتاب (mode='reflect') که پیکسلهای لبه را بازمیتاباند—معادلِ گسستهٔ همان روشِ تصویرها.
نمای گسسته: ماتریسِ توپلیتس و کد
در رایانه، سیگنال یا تصویر به یک بردار یا شبکهٔ گسسته نمونهبرداری میشود و انتگرالِ کانولوشن به یک مجموعِ گسسته بدل میشود:
این عمل را میتوان بهصورتِ یک ضربِ ماتریسی نوشت. اگر سیگنالِ ورودی بردارِ \(\vec{f} = [f_1, f_2, f_3, f_4]^\top\) و فیلترِ گاوسیِ گسسته (تابعِ گرین) دارای وزنهای \([g_{-1}, g_0, g_1]\) باشد، با فرضِ صفرگذاری در مرزها، \(u = M\vec{f}\) چنین است:
این ماتریسِ راهراهِ قطری، ماتریسِ توپلیتس (Toeplitz) نام دارد؛ ضرب در آن دقیقاً معادلِ اجرای یک حلقهٔ فیلترِ ۱بعدی است. در پردازشِ تصویرِ ۲بعدی، این به ماتریسِ بلوکیتوپلیتسِبلوکهایتوپلیتس (BTTB) گسترش مییابد.
البته در عمل، ساختنِ این ماتریسِ بزرگ (که بیشترش صفر است) ناکارآمد است؛ بهجای آن از توابعِ بهینهٔ پنجرهٔ لغزان استفاده میکنیم:
import numpy as np
from scipy.ndimage import convolve
# a 1D input signal (e.g., a sudden temperature spike)
f = np.array([0.0, 0.0, 10.0, 0.0, 0.0])
# the discrete Green's function (Gaussian blur kernel); sums to 1 to conserve "heat"
G = np.array([0.25, 0.5, 0.25])
u = convolve(f, G, mode='constant', cval=0.0)
print(u) # -> [0.0, 2.5, 5.0, 2.5, 0.0]
# the sharp spike at index 2 has diffused to its neighbors
حوزهٔ بسامد: قضیهٔ کانولوشن و میانبُرِ FFT
محاسبهٔ کانولوشن در حوزهٔ مکان (یا ضربِ ماتریسی) با بزرگشدنِ هسته بهشدت کند میشود. قضیهٔ کانولوشن این کار را بهکلی دور میزند: کانولوشن در حوزهٔ مکان/زمان، برابرِ ضربِ نقطهبهنقطه در حوزهٔ بسامد است:
برای مثالِ گاوسیِ ما، تبدیلِ فوریهٔ یک گاوسی، باز هم یک گاوسی است: \(\mathcal{F}\{e^{-x^2/2\sigma^2}\} \propto e^{-\omega^2\sigma^2/2}\). توجه کنید که با بزرگشدنِ بسامدِ \(\omega\) (نوفهٔ پربسامد، لبههای تیز)، مقدارِ \(e^{-\omega^2}\) بهسرعت به صفر میرود. پس ضربِ \(F(\omega)\cdot G(\omega)\) بسامدهای بالا را پاک و بسامدهای پایین را عبور میدهد—یعنی یک فیلترِ پایینگذر؛ و این، همان دلیلی است که پخشِ گاوسی نوفه را هموار میکند. در عمل، رایانهها مسئلهٔ فیلترِ تابعِ گرین را با میانبُرِ زیر حل میکنند:
که برای سیگنالهای بزرگ، بهمراتب سریعتر از ضرب در یک ماتریسِ توپلیتسِ غولپیکر است (همان ایدهای که در فصلِ حوزهٔ بسامد هم دیدیم).
پیوند با هوشِ مصنوعی: CNN، FNO و PINN
در فیزیکِ کلاسیک، انسانها ساعتها صرفِ حلِ معادلات میکنند تا تابعِ گرین (مانندِ گاوسی یا پتانسیلِ کولن) را بیابند. در هوشِ مصنوعی، این تابعها مستقیماً از داده آموخته میشوند.
شبکههای عصبیِ پیچشی (CNN). هر لایهٔ یک CNN فیلترهای کوچکی را روی ورودی میلغزاند—که اکنون میدانیم همان کانولوشنِ گسسته است. اگر یک CNN برای پیشبینیِ رفتارِ یک سامانهٔ فیزیکی (مثلاً جریانِ سیال یا میدانِ الکتریکی) آموزش ببیند، وزنهای فیلترهای پیچشیاش عملاً به تابعهای گرینِ آن سامانه تبدیل میشوند.
عملگرهای عصبیِ فوریه (FNO). این خانواده، رابطه را به حدِ نهایی میرساند: بهجای کار در حوزهٔ مکان، ورودی را به حوزهٔ بسامد میبرند، در مجموعهای از وزنهای آموختهشده (همان تابعِ گرینِ بسامدی) ضرب میکنند، و با FFT وارون بازمیگردانند—دقیقاً همان دستورِ \(\text{InverseFFT}(\text{FFT}(f)\times\text{FFT}(G))\)، اما با هستهای که از داده آموخته شده. به این ترتیب، شبکه یک «فیلترِ» جهانی میآموزد که معادلاتِ پیچیدهٔ مهندسی را هزاران بار سریعتر از شبیهسازهای سنتی حل میکند، چون قوانینِ فیزیک را به یک مسئلهٔ فیلترِ دادهمحور بدل میکند.
شبکههای عصبیِ آگاه از فیزیک (PINN). رویکردی مکمل: بهجای آموختنِ تابعِ گرین از داده، خودِ معادلهٔ دیفرانسیل را بهعنوانِ یک جریمه در تابعِ زیان میگنجانند. برای شبکهای که خروجیِ \(u_\theta(x,t)\) را پیشبینی میکند، پسماند (residual) را تعریف میکنند:
اگر شبکه از فیزیک پیروی کند، \(R\) صفر است؛ وگرنه جریمه میشود. تابعِ زیانِ کل سه بخش دارد: \(\mathcal{L} = \mathcal{L}_{\text{data}} + \mathcal{L}_{\text{boundary}} + \mathcal{L}_{\text{physics}}\)، که بهترتیب شبکه را به تطبیق با اندازهگیریها، رعایتِ شرایطِ مرزی، و اطاعت از معادلهٔ دیفرانسیل وادار میکنند. مشتقهای لازم با مشتقگیریِ خودکار (autodiff)—همان موتورِ پسانتشار—دقیق و تحلیلی محاسبه میشوند. با کمینهکردنِ این زیان، بهینهسازی مانندِ یک پیکرتراش، تمامِ توابعی را که فیزیک را نقض میکنند کنار میگذارد تا وزنهای شبکه به تابعِ گرینِ پیوستهٔ سامانه نزدیک شوند.
جمعبندی
در این فصل، فیلترها را ساختیم: ابزارهایی که بسامدهای خاصی را عبور میدهند و بقیه را حذف میکنند. فیلترهای پایینگذر، بالاگذر، میانگذر و میاننگذر هر کدام کاربردِ خاصِ خود را دارند، از حذفِ نوفهٔ خطِ برق تا جداکردنِ باندهای مغزی. دیدیم که هر فیلتر با پاسخِ بسامدیِ خود توصیف میشود، و چگونه با scipy.signal فیلترهای باترورث را طراحی و اعمال کنیم. و سرانجام دیدیم که فیلتر کردن در اصل همان کانولوشنِ فصلِ پیش است، که در حوزهٔ بسامد به ضربِ ساده بدل میشود. در فصلِ بعد، به تحلیلِ زمان–بسامد میپردازیم، جایی که هم زمان و هم بسامد را همزمان دنبال میکنیم.