کانولوشن با FFT در پایتون؛ صفرگذاری، طول فرد و آزمون مرزها

برابر بودن طولِ دو آرایه، برای پذیرش یک زنجیرهٔ پردازش کافی نیست. ممکن است تبدیل مستقیم، ضرب طیف و تبدیل معکوس هرکدام درست اجرا شوند، اما اتصال آن‌ها قراردادِ انتهای دنباله را عوض کند. در این درس، یک ورودیِ ولتاژ ساختگی و ضرایب بی‌بعد داریم؛ کانولوشن خطی را با حاصل‌ضرب FFT مقایسه می‌کنیم و نشان می‌دهیم صفرگذاری چه خطایی را رفع می‌کند. مثال، ثبت حسگر، آزمایش سخت‌افزار یا سنجش سرعت نیست. مسئلهٔ اصلی، آزمونِ یکپارچگی همین سه مرحله و نگه‌داشتن انتظارِ مستقل از پیاده‌سازی است.

۱. انتهای دنباله، بخشی از قرارداد است

برای دو دنبالهٔ محدود، قرارداد خطی این است که بیرون بازهٔ هر ورودی مقدار صفر فرض شود. خروجیِ اندیس k جمعِ حاصل‌ضرب‌های x[i] و h[k−i] است؛ سهمی که اندیسِ ضریبش بیرون بازه باشد وارد جمع نمی‌شود. اگر طول‌ها P و Q باشند، آخرین خانهٔ ممکن اندیس P+Q−2 دارد؛ بنابراین خروجی کامل P+Q−1 خانه می‌خواهد. مستندات کانولوشن خطی و حالت full در NumPy همین قراردادِ طول را بیان می‌کند. صفرِ بیرون بازه، فرضِ ریاضی مدل است؛ دادهٔ دریافت‌نشده را اندازه‌گیریِ صفر معرفی نمی‌کنیم.

ورودیِ تألیفی x برابر یک، دو و سه ولت، و ضرایب h برابر یک و منفی یک‌اند. محاسبهٔ دستی چهار خروجیِ یک، یک، یک و منفی سه ولت می‌دهد. خروجی دوم از دو منهای یک، و خروجی سوم از سه منهای دو ساخته می‌شود. خانهٔ چهارم هنوز لازم است، چون ضریب منفی یک در آخرین نمونه ضرب می‌شود. این دنباله با انتخاب سه خروجیِ دلخواه از میانه یا ابتدای پاسخ یکسان نیست.

در آزمون بازگشت FIR علّی، انتقالِ حافظه میان بسته‌ها و تخلیهٔ صریحِ انتهای پاسخ مطرح بود. اینجا کل دو دنباله از ابتدا در حافظه‌اند و حالتِ جریانی نداریم. همچنین معیار دامنهٔ FFT ولتاژ به تفسیر طیف یک‌طرفه مربوط است؛ دو برابر کردنِ خانه‌های طیف دامنه را نباید به ضربِ طیف‌های این الگوریتم منتقل کرد.

۲. صفرگذاری و طولِ لازم برای خروجی خطی

DFT با طول N، اندیس‌ها را در یک دورهٔ Nتایی می‌بیند. حاصلِ تبدیل معکوسِ ضرب دو DFT، کانولوشن حلقویِ همان دوره است. وقتی N حداقل برابر طولِ هر ورودی باشد، ورودی‌ها بریده نمی‌شوند؛ اما این شرط به‌تنهایی جلوی برگشتِ انتهای خروجی به ابتدای آن را نمی‌گیرد. برای مثالِ ما، N برابر سه ورودی‌ها را کامل نگه می‌دارد، ولی پاسخ خطی چهار خانه دارد. مقدار منفی سهِ خانهٔ آخر به خانهٔ صفر برمی‌گردد و نتیجه منفی دو، یک و یک می‌شود.

برای بازیابی پاسخ خطی، هر دو ورودی را تا N حداقل برابر P+Q−1 با صفر ادامه می‌دهیم. در این حالت هیچ سهمِ غیرصفرِ پاسخ کامل از انتهای دوره عبور نمی‌کند. توضیح padding برای جلوگیری از کانولوشن حلقوی مبنای این انتخاب است. اگر N بزرگ‌تر باشد، نخست P+Q−1 خروجی را برمی‌داریم؛ خانه‌های اضافی در مدلِ دقیق صفرند و در محاسبهٔ شناور ممکن است باقیماندهٔ بسیار کوچک داشته باشند. این با بریدنِ ورودیِ بلند توسط FFT فرق دارد.

مقیاس نیز بخشی از قرارداد است. طبق تعریف DFT و نرمال‌سازی backward، تبدیل مستقیم بی‌مقیاس است و معکوس ضریب یک تقسیم‌بر N دارد. کد همین گزینه را در هر سه فراخوانی صریح می‌نویسد؛ تقسیم دوبارهٔ خروجی بر N لازم نیست. نمونهٔ x واحد ولت دارد، h بی‌بعد است و جمع حاصل‌ضرب‌ها همچنان ولت می‌دهد. این جمع گسسته است، نه انتگرالِ هم‌نهشتیِ پیوسته؛ در فاصلهٔ زمانی ضرب نمی‌کنیم و بدون نرخ نمونه‌برداری، محور هرتز یا زمان فیزیکی نمی‌سازیم.

۳. کد کامل با طولِ صریح و ورودیِ محدود

برای کوتاه و روشن ماندن قرارداد، تابع فقط ndarray واقعی با dtype دقیقِ float64، یک بُعد و یک تا ۳۲ عضو می‌پذیرد. هر عضو باید متناهی و قدرمطلقش حداکثر هزار باشد؛ این حد برای هر دو آرایه اعمال می‌شود. فهرست پایتونی، آرایهٔ صحیح، bool، نوع مختلط و float32 خودکار تبدیل نمی‌شوند. فراخواننده باید تبدیلِ موردنظرش را پیش از ورود صریح انجام دهد. این محدودیتِ آموزشیِ برنامهٔ ماست، نه محدودیت عمومی NumPy.

طول تبدیل باید int دقیق پایتونی و حداکثر ۱۲۸ باشد؛ bool و اسکالر صحیح NumPy نیز در این رابط رد می‌شوند. تابع حلقوی N کوچک‌تر از طولِ ورودی‌ها را نمی‌پذیرد، چون rfft در آن حالت ورودی را می‌بُرد. تابع خطی شرطِ قوی‌ترِ N حداقل برابر طول پاسخ کامل را اعمال می‌کند. وقتی N داده نشود، همان طول کامل انتخاب می‌شود؛ لازم نیست طولِ مثال را به توانِ دو گرد کنیم.

پارامتر n در irfft عمداً حذف نشده است. قرارداد طول خروجیِ irfft می‌گوید اگر طیف m عضو داشته باشد، طولِ پیش‌فرض ۲(m−۱) است. طیفِ حاصل از پنج نمونه، سه عضو دارد؛ همان اندازه می‌تواند متعلق به خروجیِ چهارعضوی هم باشد. پس اندازهٔ طیف، زوج یا فرد بودنِ طول اصلی را به‌تنهایی تعیین نمی‌کند. ورودی مختلط نیز پیش از rfft رد می‌شود تا حذفِ بی‌صدای بخش موهومی رخ ندهد.

فایل زیر با Python 3.14.5 و NumPy 2.5.4 اجرا شده است؛ صفحات جاری، راهنمای NumPy v2.5 را نمایش می‌دهند. آرایه‌های ورودی کپی می‌شوند و محاسبه فایل یا شبکه ندارد. خروجیِ JSON، نسخه‌ها و بررسی‌ها را همراه نمونه نگه می‌دارد.

"""Bounded offline FFT convolution: finite synthetic volt samples, no benchmark."""
import json
import platform
import numpy as np

MAX_SAMPLES = 32
MAX_FFT = 128
MAX_ABS = 1000.0
ATOL = 1e-6  # Own test tolerance in output units, not sensor accuracy.


def vector(value):
    """Own API: exact ndarray, native float64, finite 1D, 1..32 samples."""
    if type(value) is not np.ndarray or value.dtype != np.dtype('float64'):
        raise TypeError('Expected an ndarray with native float64 dtype')
    if value.ndim != 1 or not 1 <= value.size <= MAX_SAMPLES:
        raise ValueError('Expected one dimension and 1..32 samples')
    if not np.isfinite(value).all() or np.any(np.abs(value) > MAX_ABS):
        raise ValueError('Samples must be finite and abs <=1000')
    return value.copy()


def length(n, minimum):
    if type(n) is not int:  # bool and numpy integer scalars are not this API.
        raise TypeError('n must be an exact Python int')
    if not minimum <= n <= MAX_FFT:
        raise ValueError('n must preserve both inputs and fit the teaching cap')
    return n


def spectral_product(a, b, n):
    first = np.fft.rfft(a, n=n, norm='backward')
    second = np.fft.rfft(b, n=n, norm='backward')
    return np.fft.irfft(first * second, n=n, norm='backward')


def circular_fft(x, h, n):
    """Return n-period circular convolution without silently cropping inputs."""
    a, b = vector(x), vector(h)
    return spectral_product(a, b, length(n, max(a.size, b.size)))


def linear_fft(x, h, n=None):
    """Return full linear convolution, length len(x)+len(h)-1."""
    a, b = vector(x), vector(h)
    full_length = a.size + b.size - 1
    n = full_length if n is None else length(n, full_length)
    return spectral_product(a, b, n)[:full_length]


def direct_linear(x, h):
    """Independent time-domain scatter sum: no FFT or numpy.convolve."""
    answer = [0.0] * (len(x) + len(h) - 1)
    for i, first in enumerate(x):
        for j, second in enumerate(h):
            answer[i + j] += first * second
    return answer


def folded_linear(x, h, n):
    """Independent reference: fold the direct linear result modulo n."""
    answer = [0.0] * n
    for index, value in enumerate(direct_linear(x, h)):
        answer[index % n] += value
    return answer


def close(actual, expected):
    actual, expected = np.asarray(actual), np.asarray(expected)
    return actual.shape == expected.shape and bool(
        np.allclose(actual, expected, rtol=0.0, atol=ATOL))


def rejected(kind, operation):
    try:
        operation()
    except kind:
        return True
    return False


def main():
    checks = {}

    def check(name, passed):
        assert passed, name
        checks[name] = True

    x = np.array([1.0, 2.0, 3.0], dtype=np.float64)
    h = np.array([1.0, -1.0], dtype=np.float64)
    check('hand_linear_four_samples', close(linear_fft(x, h), [1, 1, 1, -3]))
    check('hand_circular_wrap_n3', close(circular_fft(x, h, 3), [-2, 1, 1]))

    linear_errors, circular_errors, library_agreements = [], [], []
    for nx in range(1, 9):
        for nh in range(1, 9):
            # Quarter-integer values are exactly representable in binary.
            a = np.array([((j + 2) * (nx + 1) % 11 - 5) / 4
                          for j in range(nx)], dtype=np.float64)
            b = np.array([((j + 3) * (nh + 2) % 9 - 4) / 4
                          for j in range(nh)], dtype=np.float64)
            expected = direct_linear(a.tolist(), b.tolist())
            linear_errors.append(float(np.max(np.abs(linear_fft(a, b) - expected))))
            n = max(nx, nh)
            folded = folded_linear(a.tolist(), b.tolist(), n)
            circular_errors.append(float(np.max(np.abs(circular_fft(a, b, n) - folded))))
            library_agreements.append(close(linear_fft(a, b), np.convolve(a, b, mode='full')))
    check('direct_time_domain_oracle_64_cases', max(linear_errors) <= ATOL)
    check('explicit_modulo_folding_oracle_64_cases', max(circular_errors) <= ATOL)
    check('numpy_full_convolve_secondary_comparison', all(library_agreements))

    odd_h = np.array([1.0, 0.0, -1.0], dtype=np.float64)
    product = np.fft.rfft(x, n=5) * np.fft.rfft(odd_h, n=5)
    omitted = np.fft.irfft(product)  # Deliberate counterexample, defaults to4.
    explicit = np.fft.irfft(product, n=5)
    check('odd_irfft_explicit_length_preserves_five',
          len(omitted) == 4 and close(explicit, [1, 2, 2, -2, -3]))
    check('extra_zero_padding_keeps_linear_prefix_and_zero_tail',
          close(circular_fft(x, h, 7), [1, 1, 1, -3, 0, 0, 0])
          and close(linear_fft(x, h, 7), [1, 1, 1, -3]))
    unit = np.array([1.0], dtype=np.float64)
    zero = np.zeros(3, dtype=np.float64)
    check('impulse_zero_and_gain_have_expected_scale',
          close(linear_fft(x, unit), x)
          and close(linear_fft(x, zero), [0, 0, 0, 0, 0])
          and close(linear_fft(2 * x, h), [2, 2, 2, -6]))
    largest = np.full(32, MAX_ABS, dtype=np.float64)
    check('maximum_input_length_magnitude_and_fft128',
          close(linear_fft(largest, largest, 128), direct_linear(largest.tolist(), largest.tolist())))
    before_x, before_h = x.copy(), h.copy()
    circular_fft(x, h, 3)
    linear_fft(x, h, 7)
    check('caller_arrays_not_modified', np.array_equal(x, before_x) and np.array_equal(h, before_h))

    bad_shapes = [np.array([], dtype=np.float64), np.zeros((1, 2), dtype=np.float64),
                  np.array(1.0, dtype=np.float64), np.zeros(33, dtype=np.float64)]
    check('empty_scalar_multidimensional_and_length33_rejected',
          all(rejected(ValueError, lambda value=value: linear_fft(value, h)) for value in bad_shapes))
    bad_types = [[1.0], np.array([1], dtype=np.int64), np.array([True]),
                 np.array([1j]), np.array(['1']), np.array([1], dtype=object),
                 np.array([1], dtype=np.float32)]
    check('container_integer_bool_complex_string_object_float32_rejected',
          all(rejected(TypeError, lambda value=value: linear_fft(value, h)) for value in bad_types))
    bad_values = [np.array([value], dtype=np.float64)
                  for value in (np.nan, np.inf, -np.inf, 1000.1, -1000.1)]
    check('nonfinite_and_out_of_bound_samples_rejected',
          all(rejected(ValueError, lambda value=value: linear_fft(value, h)) for value in bad_values))
    bad_n_types = [True, 4.0, np.int64(4), None]
    check('transform_length_type_minimum_and_teaching_cap_rejected',
          all(rejected(TypeError, lambda n=n: circular_fft(x, h, n)) for n in bad_n_types)
          and all(rejected(ValueError, lambda n=n: circular_fft(x, h, n)) for n in (-1, 0, 2, 129))
          and rejected(ValueError, lambda: linear_fft(x, h, 3)))
    assert len(checks) == 14
    report = {
        'python': platform.python_version(), 'numpy': np.__version__, 'checks': checks,
        'meaningful_named_groups': len(checks), 'family_cases_per_oracle': 64,
        'example': {'x_volts': x.tolist(), 'h_dimensionless': h.tolist(),
                    'linear_volts': np.round(linear_fft(x, h), 10).tolist(),
                    'circular_n3_volts': np.round(circular_fft(x, h, 3), 10).tolist(),
                    'linear_length': 4, 'odd_explicit_length': len(explicit),
                    'odd_omitted_length': len(omitted)},
        'observed_max_abs_error_direct_family': max(linear_errors),
        'observed_max_abs_error_folded_family': max(circular_errors),
        'absolute_test_tolerance_volts': ATOL,
        'scope': 'Synthetic discrete volt samples and dimensionless coefficients; no timing, sensor, continuous integral, stream-state or universal accuracy claim.'}
    print(json.dumps(report, ensure_ascii=False, indent=2))


if __name__ == '__main__':
    main()

۴. شاهدِ دستی، مرجعِ مستقیم و آزمونِ حلقوی

مرجعِ اصلی در حوزهٔ اندیس کار می‌کند: دو حلقهٔ مستقل، سهم هر جفت ورودی و ضریب را در خانهٔ مجموع اندیس‌ها می‌ریزند. این مسیر از FFT، تبدیل معکوس و convolve کتابخانه استفاده نمی‌کند. برای شاهد حلقوی، همان پاسخ مستقیم را با باقیماندهٔ تقسیمِ اندیس بر N تا می‌کنیم. این مرجع، برگشتِ انتهای پاسخ را آشکار می‌کند؛ مقایسهٔ تبدیل و معکوسِ یکدیگر به‌تنهایی ممکن است دو خطای مقیاسِ سازگار را پنهان کند.

برای هر جفت طول از یک تا هشت، داده‌های کوچکِ ربع‌صحیح می‌سازیم؛ در مجموع ۶۴ حالت برای مرجع خطی و همان تعداد برای مرجع حلقوی داریم. این‌ها دو گروه بررسی‌اند، نه ۱۲۸ مقاله یا گروه مستقل. مقایسه با numpy.convolve در حالت full، کنترلِ سوم و مکمل است. آزمونِ پاسخ ضربه، ورودی صفر، تغییر بهره، بیشینهٔ طول و دامنهٔ مجاز، دست‌نخوردن ورودی و خطاهای قرارداد نیز جدا ثبت شده‌اند.

هر چهارده گروه موفق شد. بیشترین اختلاف مشاهده‌شده در خانوادهٔ خطی حدود ۸٫۸۸ ضرب‌در ده به توان منفی شانزده ولت و در خانوادهٔ حلقوی حدود ۴٫۴۴ ضرب‌در ده به توان منفی شانزده ولت بود. تلورانس مطلق آزمون یک میکروولت است و خطای نسبی در مقایسه صفر انتخاب شده؛ این آستانه دقت حسگر نیست و از این داده‌های محدود، کران عمومیِ خطای FFT نتیجه نمی‌گیریم. در شاهدِ طول فرد، حذف n چهار خروجی می‌دهد، درحالی‌که طول صریح پنج، پاسخِ پنج‌عضوی مورد انتظار را بازمی‌گرداند.

{
  "python": "3.14.5",
  "numpy": "2.5.4",
  "checks": {
    "hand_linear_four_samples": true,
    "hand_circular_wrap_n3": true,
    "direct_time_domain_oracle_64_cases": true,
    "explicit_modulo_folding_oracle_64_cases": true,
    "numpy_full_convolve_secondary_comparison": true,
    "odd_irfft_explicit_length_preserves_five": true,
    "extra_zero_padding_keeps_linear_prefix_and_zero_tail": true,
    "impulse_zero_and_gain_have_expected_scale": true,
    "maximum_input_length_magnitude_and_fft128": true,
    "caller_arrays_not_modified": true,
    "empty_scalar_multidimensional_and_length33_rejected": true,
    "container_integer_bool_complex_string_object_float32_rejected": true,
    "nonfinite_and_out_of_bound_samples_rejected": true,
    "transform_length_type_minimum_and_teaching_cap_rejected": true
  },
  "meaningful_named_groups": 14,
  "family_cases_per_oracle": 64,
  "example": {
    "x_volts": [
      1.0,
      2.0,
      3.0
    ],
    "h_dimensionless": [
      1.0,
      -1.0
    ],
    "linear_volts": [
      1.0,
      1.0,
      1.0,
      -3.0
    ],
    "circular_n3_volts": [
      -2.0,
      1.0,
      1.0
    ],
    "linear_length": 4,
    "odd_explicit_length": 5,
    "odd_omitted_length": 4
  },
  "observed_max_abs_error_direct_family": 8.881784197001252e-16,
  "observed_max_abs_error_folded_family": 4.440892098500626e-16,
  "absolute_test_tolerance_volts": 1e-06,
  "scope": "Synthetic discrete volt samples and dimensionless coefficients; no timing, sensor, continuous integral, stream-state or universal accuracy claim."
}

۵. چهار پرسش دربارهٔ طول و مقیاس

پرسش اول: آیا پاسخ حلقویِ سه‌عضوی فقط سه خانهٔ اولِ پاسخ خطی است؟ پاسخ: خیر؛ منفی سهِ انتهای پاسخ خطی به خانهٔ صفر افزوده می‌شود و یک به منفی دو تغییر می‌کند. این خطا با حذفِ خانهٔ آخر از نمایش رفع نمی‌شود؛ باید قراردادِ دوره و طول تبدیل درست انتخاب شود.

پرسش دوم: برای طول‌های سه و سه، کمترین N مجاز در تابع خطی چیست؟ پاسخ: پنج. کد، طول پنج را پشتیبانی می‌کند و همان پنج را به irfft می‌دهد. استفاده از پیش‌فرضِ معکوس، با وجود اندازهٔ درستِ طیف، طول چهار می‌سازد. طول هشت هم مجاز است، ولی الزامِ ریاضیِ این مثال نیست.

پرسش سوم: اگر همهٔ نمونه‌های ورودی دو برابر شوند، آیا باید خروجی را دوباره بر N تقسیم کنیم؟ پاسخ: خیر؛ با ضرایب ثابت، خروجی دو برابر می‌شود و شاهدِ ما دو، دو، دو و منفی شش ولت خواهد بود. ضریبِ معکوس در قرارداد backward از قبل اعمال شده است.

پرسش چهارم: آیا عبور از این آزمون‌ها، دریافتِ پیوسته یا سریع‌تر بودن روش را ثابت می‌کند؟ پاسخ: خیر؛ این برنامه حالتِ بسته‌ایِ محدود دارد و زمان اجرا اندازه‌گیری نشده است. برای پردازش جریانی باید سیاستِ حافظه و مرز بلوک جدا تعریف و آزموده شود. منابع در ۲۷ مهر ۱۴۰۵ خوانده شدند؛ تصاویرِ همراه، مفهومی‌اند و نمایشِ داده یا نتیجهٔ واقعی آزمایش نیستند.

منابع و مطالعهٔ بیشتر

• NumPy v2.5؛ قرارداد DFT و نرمال‌سازی backward
• NumPy v2.5؛ کانولوشن خطی، طول full و صفرگذاری
• NumPy v2.5؛ طول صریح irfft و تفاوتِ زوج و فرد
• NumPy v2.5؛ طول rfft، صفرگذاری و بریدن ورودی

بررسی منابع: ۲۷ مهر ۱۴۰۵. تحلیل‌ها و پیشنهادهای اجرایی، جمع‌بندی نویسنده‌اند.