بازسازی STFT با پایتون؛ پنجره، گام و وزن هم‌پوشانی صوت ساختگی

در تحلیل صوتِ متغیر با زمان، طیفِ یکپارچهٔ کل قطعه، نمایش زمان‌محلیِ تغییرات را به‌تنهایی در اختیارمان نمی‌گذارد. می‌توان صوت را به فریم‌های کوتاه تقسیم و هر فریم را تبدیل کرد؛ اما اگر بعداً بخواهیم نمونه‌ها را بازسازی کنیم، جمع سادهٔ فریم‌ها کافی نیست. این درس تبدیل زمان‌کوتاه و بازسازی وزن‌دار را روی یک سیگنال کوچکِ ساختگی اجرا می‌کند تا نقش پنجره، گام و پوشش مرز روشن شود.

۱. صوت مصنوعی را به فریم‌های پنجره‌خورده تقسیم کنیم

ورودی هشتاد نمونهٔ حقیقی با نرخ ۸۰۰۰ نمونه در ثانیه است؛ مدت کل ده میلی‌ثانیه. چهل نمونهٔ اول سینوس ۵۰۰ هرتز با دامنهٔ ۰٫۴ و بقیه سینوس ۱۵۰۰ هرتز با دامنهٔ ۰٫۲۵ دارند؛ به هر دو آفست ۰٫۰۵ اضافه شده است. دامنه بی‌واحد است و این قطعه نه گفتار طبیعی است، نه ولتاژ حسگر یا نمونهٔ ضبط تلفنی.

طول فریم شانزده و گام هشت نمونه انتخاب شده‌اند؛ یعنی فریم دو میلی‌ثانیه و فاصلهٔ شروع‌ها یک میلی‌ثانیه است. هر فریم در پنجرهٔ متقارن Hann ضرب می‌شود. فرمول hanning در NumPy برای طول شانزده مخرج پانزده دارد و دو سرش صفر است؛ آن را با نوع تناوبیِ درس قبلی یکی نمی‌گیریم.

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

۲. هم‌پوشانی به وزن بازسازی نیاز دارد

تبدیل جلو روی آخرین محورِ آرایهٔ فریم‌ها انجام می‌شود. در rfft، برای ورودی حقیقی با طول شانزده، نُه ضریبِ فرکانسی کافی است؛ بخش منفی از تقارن به دست می‌آید. در برنامهٔ ما محور اول زمانِ فریم و محور دوم فرکانس است، پس شکل خروجی یازده در نُه خواهد بود.

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

توضیح شرط NOLA در SciPy غیرصفر بودن همین مخرج را بیان می‌کند. این منبع برای رابطهٔ نظری خوانده شده است؛ SciPy در محیط فعلی نصب نیست و اجرا نشده است. فریم‌بندی و بازسازی حاضر، کد تألیفیِ مبتنی بر NumPy است و پیاده‌سازیِ API خاص SciPy معرفی نمی‌شود.

۳. پوشش مرز و صفر مخرج را جدا بررسی کنیم

پیش و پس از سیگنال، هشت نمونهٔ صفر قرار می‌دهیم تا نمونه‌های مرزی نیز زیر قسمت غیرصفر پنجره قرار بگیرند. سپس بعد از بازسازی فقط هشتاد نمونهٔ اصلی را برمی‌داریم. صفرها دادهٔ واقعیِ قبل و بعد تماس نیستند؛ قرارداد مرزیِ آزمایش‌اند. برای طول‌های دیگر، انتهای آرایه به اندازهٔ لازم برای فریم کامل هم پر می‌شود.

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

دو میلی‌ثانیه بودن فریم، تأخیر کل یک سامانه را تعیین نمی‌کند. این برنامه قطعهٔ کامل را از قبل دارد؛ انتظار برای نمونهٔ آینده، زمان محاسبه و انتقال در آن مدل نشده‌اند. اجرای آفلاین را نباید پردازش زنده با تضمین موعد نامید.

Hann بدون هم‌پوشانی نیز حفرهٔ وزن دارد، ولی پنجرهٔ مستطیلیِ بدون هم‌پوشانی در همین آزمایش بازسازی می‌شود. پس حکم «همیشه هم‌پوشانی لازم است» دقیق نیست؛ ترکیب پنجره، گام و پوشش مهم است. برای پرهیز از تقسیم عددیِ نامطمئن، در قرارداد ما وزنِ درخواستی باید از 1e-14 بزرگ‌تر باشد؛ این آستانه با تلورانسِ تابع SciPy یکی نیست.

۴. کد را با DFT مستقیم و پاسخ دستی بسنجیم

فایل stft_overlap.py با Python 3.14.5 و NumPy 2.4.6 اجرا شده است. حالت نرمال‌سازی backward صریح است: تبدیل جلو مقیاس نمی‌شود و تبدیل معکوس عامل طول را لحاظ می‌کند. در irfft، طول شانزده را صریح می‌دهیم تا تفسیر طول و محور پنهان نباشد.

"""Original NumPy STFT/weighted overlap-add fixture; SciPy is not executed."""
import cmath
import math
import json
import platform
import numpy as np

ATOL = 2e-13
MIN_WEIGHT = 1e-14


def real_vector(value, minimum, maximum):
    if type(value) is not np.ndarray or value.dtype != np.float64 or value.ndim != 1:
        raise ValueError('one-dimensional actual float64 ndarray required')
    if not minimum <= len(value) <= maximum or not np.isfinite(value).all():
        raise ValueError('finite bounded vector required')


def analyse(signal, window, hop, *, pad_edges=True):
    real_vector(signal, 1, 256)
    real_vector(window, 4, 32)
    length = len(window)
    if type(hop) is not int or not 1 <= hop <= length or type(pad_edges) is not bool:
        raise ValueError('invalid hop or padding policy')
    if np.any(window < 0):
        raise ValueError('own profile requires nonnegative analysis weights')
    left = length - hop if pad_edges else 0
    right = length - hop if pad_edges else 0
    padded = np.pad(signal, (left, right))
    extra = max(0, length - len(padded))
    extra += (-(len(padded) + extra - length)) % hop
    padded = np.pad(padded, (0, extra))
    starts = list(range(0, len(padded) - length + 1, hop))
    frames = np.stack([padded[start:start + length] * window for start in starts])
    spectrum = np.fft.rfft(frames, n=length, axis=-1, norm='backward')
    return {'spectrum': spectrum, 'frames': frames, 'starts': starts,
            'window': window.copy(), 'hop': hop, 'signal_length': len(signal),
            'left_padding': left, 'right_padding': right + extra,
            'padded_length': len(padded)}


def reconstruct(plan, spectrum=None):
    z = plan['spectrum'] if spectrum is None else spectrum
    if type(z) is not np.ndarray or z.dtype != np.complex128 or z.shape != plan['spectrum'].shape or not np.isfinite(z).all():
        raise ValueError('finite complex128 spectrum with original shape required')
    length = len(plan['window'])
    frames = np.fft.irfft(z, n=length, axis=-1, norm='backward')
    numerator = np.zeros(plan['padded_length'], dtype=np.float64)
    denominator = np.zeros_like(numerator)
    for row, start in enumerate(plan['starts']):
        numerator[start:start + length] += frames[row] * plan['window']
        denominator[start:start + length] += plan['window'] ** 2
    left = plan['left_padding'];right = left + plan['signal_length']
    weights = denominator[left:right]
    if np.any(weights <= MIN_WEIGHT):
        raise ValueError('uncovered or numerically unsafe requested reconstruction samples')
    return numerator[left:right] / weights, weights.copy()


def direct_dft(frame):
    # Independent finite sum; no NumPy FFT/ifft or SciPy transform.
    length = len(frame)
    return [sum(float(value) * cmath.exp(-2j * math.pi * k * n / length)
                for n, value in enumerate(frame)) for k in range(length // 2 + 1)]


def close(a, b):
    return bool(np.allclose(a, b, rtol=0, atol=ATOL))


def rejects(fn):
    try:
        fn()
    except ValueError:
        return True
    return False


def run_example():
    fs = 8000
    n = np.arange(80, dtype=np.float64)
    signal = 0.05 + np.where(n < 40, 0.4 * np.sin(2 * np.pi * 500 * n / fs),
                            0.25 * np.sin(2 * np.pi * 1500 * n / fs))
    window = np.hanning(16).astype(np.float64)
    signal_before = signal.copy();window_before = window.copy()
    plan = analyse(signal, window, 8)
    restored, weights = reconstruct(plan)
    scalar_spectra = np.array([direct_dft(frame) for frame in plan['frames']], dtype=np.complex128)
    manual_window = np.array([0.5 - 0.5 * math.cos(2 * math.pi * j / 15) for j in range(16)])
    holes = np.array([1.0, 0.0] * 8, dtype=np.float64)
    rect = analyse(signal, np.ones(16, dtype=np.float64), 16)
    short = np.array([0.3], dtype=np.float64)
    odd = np.linspace(-0.2, 0.3, 65, dtype=np.float64)
    phase_removed = np.abs(plan['spectrum']).astype(np.complex128)
    modified, _ = reconstruct(plan, phase_removed)
    zero = np.zeros(19, dtype=np.float64)
    checks = {
        'unchanged_stft_reconstructs_original': close(restored, signal),
        'all_frame_coefficients_match_direct_finite_dft': close(plan['spectrum'], scalar_spectra),
        'four_point_hand_dft_and_inverse_normalization': close(np.fft.rfft(np.array([1., 0., -1., 0.]), norm='backward'), [0, 2, 0]) and close(np.fft.irfft(np.array([0, 2, 0], dtype=np.complex128), n=4, norm='backward'), [1, 0, -1, 0]),
        'symmetric_hann_matches_scalar_formula': close(window, manual_window) and bool(window[0] == window[-1] == 0),
        'requested_boundaries_have_nonzero_weights': bool(np.min(weights) > MIN_WEIGHT) and close(restored[[0, -1]], signal[[0, -1]]),
        'no_edge_padding_exposes_first_sample_hole': rejects(lambda: reconstruct(analyse(signal, window, 8, pad_edges=False))),
        'hann_without_overlap_has_weight_holes': rejects(lambda: reconstruct(analyse(signal, window, 16))),
        'alternating_zero_window_exposes_interior_holes': rejects(lambda: reconstruct(analyse(signal, holes, 8))),
        'rectangular_nonoverlap_can_reconstruct': close(reconstruct(rect)[0], signal),
        'single_sample_and_nonmultiple_length_crop': close(reconstruct(analyse(short, window, 8))[0], short) and close(reconstruct(analyse(odd, window, 8))[0], odd),
        'zero_signal_and_linear_scaling': close(reconstruct(analyse(zero, window, 8))[0], zero) and close(reconstruct(analyse(2 * signal, window, 8))[0], 2 * restored),
        'phase_removal_is_not_original_reconstruction': bool(np.max(np.abs(modified - signal)) > 0.1),
        'spectrum_axes_and_input_preservation': plan['spectrum'].shape == (11, 9) and np.array_equal(signal, signal_before) and np.array_equal(window, window_before),
        'invalid_input_types_shapes_and_finite_values': all(rejects(lambda x=x: analyse(x, window, 8)) for x in [[], signal.astype(np.float32), signal.reshape(10, 8), np.array([np.nan]), np.array([np.inf]), np.array([1+1j])]),
        'hop_window_and_spectrum_contract': all(rejects(fn) for fn in [lambda: analyse(signal, window, True), lambda: analyse(signal, window, 0), lambda: analyse(signal, window, 17), lambda: analyse(signal, -window, 8), lambda: reconstruct(plan, plan['spectrum'][:, :-1]), lambda: reconstruct(plan, np.full((11, 9), np.nan + 0j))]),
    }
    assert all(type(x) is bool and x for x in checks.values()), checks
    return {'python': platform.python_version(), 'numpy': np.__version__,
            'dependencies': 'NumPy plus Python standard library; SciPy not installed/executed',
            'example_kind': 'Original synthetic dimensionless signal, NumPy STFT and weighted overlap-add; no speech/codec/noise reduction benchmark',
            'sample_rate_hz': fs, 'samples': len(signal), 'signal_duration_ms': len(signal) * 1000 / fs,
            'frame_length_samples': len(window), 'hop_samples': plan['hop'],
            'frame_duration_ms': len(window) * 1000 / fs, 'hop_ms': plan['hop'] * 1000 / fs,
            'frequency_bin_spacing_hz': fs / len(window), 'spectrum_shape_time_frequency': list(plan['spectrum'].shape),
            'left_padding_samples': plan['left_padding'], 'right_padding_samples': plan['right_padding'],
            'minimum_requested_weight': float(np.min(weights)),
            'maximum_reconstruction_absolute_error': float(np.max(np.abs(restored - signal))),
            'maximum_direct_dft_absolute_error': float(np.max(np.abs(scalar_spectra - plan['spectrum']))),
            'absolute_acceptance_tolerance': ATOL, 'relative_acceptance_tolerance': 0,
            'float64_epsilon': float(np.finfo(np.float64).eps),
            'phase_removed_maximum_absolute_difference': float(np.max(np.abs(modified - signal))),
            'direct_reference_coefficients_inside_one_group': int(scalar_spectra.size), 'checks': checks}


if __name__ == '__main__':
    print(json.dumps(run_example(), ensure_ascii=False, indent=2))

مرجع مستقل، برای تک‌تک ضرایب، جمع محدودِ نمایی مختلط را با cmath حساب می‌کند؛ از FFT استفاده نمی‌کند. نودونه ضریب در یک گروه مقایسه می‌شوند. یک جواب دستیِ چهارنقطه‌ای هم داریم: ورودی [1,0,-1,0] در تبدیل یک‌طرفه به [0,2,0] می‌رسد و معکوس همان چهار نمونه را برمی‌گرداند.

در اجرای ثبت‌شده، کمترین مخرجِ نمونه‌های اصلی حدود ۰٫۴۲۴۳۶ و بیشترین خطای بازسازی حدود 1.11e-16 است. اختلاف ضرایب با مرجع مستقیم حدود 7.52e-16 شد. تلورانس پذیرشِ این مثال کوچک، مطلقِ 2e-13 و نسبیِ صفر انتخاب شده؛ این حاشیه برای محاسبات محدودِ double است و کرانِ نظریِ خطای همهٔ ورودی‌ها نیست.

۵. بازسازی عددی را کیفیت صوت یا حذف نویز ننامیم

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

{
  "python": "3.14.5",
  "numpy": "2.4.6",
  "dependencies": "NumPy plus Python standard library; SciPy not installed/executed",
  "example_kind": "Original synthetic dimensionless signal, NumPy STFT and weighted overlap-add; no speech/codec/noise reduction benchmark",
  "sample_rate_hz": 8000,
  "samples": 80,
  "signal_duration_ms": 10.0,
  "frame_length_samples": 16,
  "hop_samples": 8,
  "frame_duration_ms": 2.0,
  "hop_ms": 1.0,
  "frequency_bin_spacing_hz": 500.0,
  "spectrum_shape_time_frequency": [
    11,
    9
  ],
  "left_padding_samples": 8,
  "right_padding_samples": 8,
  "minimum_requested_weight": 0.42436016005775895,
  "maximum_reconstruction_absolute_error": 1.1102230246251565e-16,
  "maximum_direct_dft_absolute_error": 7.519225344821328e-16,
  "absolute_acceptance_tolerance": 2e-13,
  "relative_acceptance_tolerance": 0,
  "float64_epsilon": 2.220446049250313e-16,
  "phase_removed_maximum_absolute_difference": 0.6089108483063609,
  "direct_reference_coefficients_inside_one_group": 99,
  "checks": {
    "unchanged_stft_reconstructs_original": true,
    "all_frame_coefficients_match_direct_finite_dft": true,
    "four_point_hand_dft_and_inverse_normalization": true,
    "symmetric_hann_matches_scalar_formula": true,
    "requested_boundaries_have_nonzero_weights": true,
    "no_edge_padding_exposes_first_sample_hole": true,
    "hann_without_overlap_has_weight_holes": true,
    "alternating_zero_window_exposes_interior_holes": true,
    "rectangular_nonoverlap_can_reconstruct": true,
    "single_sample_and_nonmultiple_length_crop": true,
    "zero_signal_and_linear_scaling": true,
    "phase_removal_is_not_original_reconstruction": true,
    "spectrum_axes_and_input_preservation": true,
    "invalid_input_types_shapes_and_finite_values": true,
    "hop_window_and_spectrum_contract": true
  }
}

اگر فقط اندازهٔ ضرایب را نگه داریم و فاز را دور بریزیم، بازسازی این نمونه با اصل فرق زیادی دارد؛ پس نمودار اندازهٔ طیف، تمام اطلاعات بازسازی نیست. هیچ حذف نویز یا بهبود تشخیص گفتار انجام نداده‌ایم. فاصلهٔ نقاط طیف ۵۰۰ هرتز است، ولی این عدد به‌تنهایی توان جداسازی دو مؤلفهٔ نزدیک را مشخص نمی‌کند.

پرسش اول: چرا بعد از جمع بر وزن تقسیم می‌کنیم؟ پاسخ: هر نمونه چند بار با وزن تحلیل و وزن بازسازی وارد شده است؛ مخرج اثر مجموع مربع این وزن‌ها را حذف می‌کند.

پرسش دوم: با مخرج صفر تقسیم را رد کنیم یا نمونه را صفر بگذاریم؟ پاسخ: این قرارداد رد می‌کند؛ دادهٔ پوشش‌نداده را نباید بازسازی‌شده نامید.

پرسش سوم: یازده فریم یعنی مدت صوت بیشتر شده؟ پاسخ: خیر؛ هم‌پوشانی و قرارداد مرز فریم‌های بیشتری می‌سازند، اما خروجی به همان هشتاد نمونه و ده میلی‌ثانیه بریده می‌شود.

پرسش چهارم: این روش همان تبدیل نرخ نمونه‌برداری است؟ پاسخ: خیر؛ نرخ در این درس ثابت مانده است. منابع در ۲۵ مهر ۱۴۰۵ خوانده شدند؛ نسخه‌های مستندات و اجرای واقعی جدا ثبت‌اند و نتیجهٔ صوت واقعی یا خبر تازه ادعا نشده است.

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

• NumPy2.4؛ rfft و محور و نرمال‌سازی
• NumPy2.4؛ irfft و طول بازسازی
• NumPy2.4؛ پنجرهٔ متقارن hanning
• SciPy؛ شرط NOLA و مجموع مربع پنجره‌ها

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