2025, Oct 18 02:32

DFT में ज़ीरो पैडिंग से खिसका DC बिन: इंडेक्सिंग व फ़्रीक्वेंसी स्केल को सही करने की गाइड

डिस्क्रीट‑टाइम सिग्नल में DFT ज़ीरो पैडिंग करते समय इंडेक्सिंग और फ़्रीक्वेंसी स्केल की गलती से DC बिन खिसकता है। कारण, सही सूत्रीकरण से समाधान जानें, विस्तार से.

डिस्क्रीट‑टाइम सिग्नल को फैलाने के लिए स्पेक्ट्रम में ज़ीरो पैडिंग करना सीधा सा लगता है, लेकिन इंडेक्स की गणित गड़बड़ा दे तो बात बिगड़ जाती है। यह संक्षिप्त मार्गदर्शिका बताती है कि फ़ूरिएर श्रेणी के गुणांकों के बीच शून्य डालने पर वास्तविक rect सिग्नल की पाँच साफ़ प्रतियों के बजाय एक कॉम्प्लेक्स, खिसका हुआ वेवफ़ॉर्म क्यों मिला—और इंडेक्सिंग व फ़्रीक्वेंसी स्केल ठीक करके इसे कैसे सुधारा जाए।

समस्या का सेटअप

काम एक साधारण rect सिग्नल a[n] से शुरू होता है, डोमेन [-1000, 1000] पर, जहाँ |n| < 100 के लिए मान 1 है और बाकी जगह 0। इसके बाद इसकी डिस्क्रीट फ़ूरिएर सीरीज़ के गुणांक निकाले जाते हैं। फिर स्पेक्ट्रम को इस तरह “फैलाया” जाता है कि हर गुणांक के बाद चार शून्य डाले जाएँ और 0.2 से स्केल किया जाए, यानी 0.2·[a0, 0, 0, 0, 0, a1, 0, 0, 0, 0, a2, …]। समय क्षेत्र में उम्मीद यह है कि मूल सिग्नल की पाँच एक जैसी प्रतियाँ मिलें, जो एक‑दूसरे से दूरी पर रखी हों।

समस्या दिखाने वाला कोड

नीचे दिया गया स्निपेट वही व्यवहार दिखाता है जहाँ इनवर्स ट्रांसफॉर्म कॉम्प्लेक्स वेवफ़ॉर्म लौटाता है और प्रतियाँ खिसकी हुई दिखती हैं। नाम समझाने के लिए हैं, लेकिन गणना की तर्कशृंखला ऊपर बताए गए तरीके का ही अनुसरण करती है।

import numpy as np
import cmath
L = 1000
JJ = complex(0, 1)
PI = np.pi
M = 2 * L + 1
# Rect सिग्नल: |n| < 100 पर मान 1
rect_sig = np.zeros(M)
for t in range(-99, 100):
    rect_sig[t + L] = 1
# छोटे FP शोर और लगभग‑पूर्णांकों को साफ़ करना
eps_cut = 1e-10
def snap_near_integers(arr, tol=eps_cut):
    r = np.real(arr)
    im = np.imag(arr)
    r[np.abs(r) < tol] = 0
    im[np.abs(im) < tol] = 0
    r_fraction = r - np.round(r)
    im_fraction = im - np.round(im)
    r[np.abs(r_fraction) < tol] = np.round(r[np.abs(r_fraction) < tol])
    im[np.abs(im_fraction) < tol] = np.round(im[np.abs(im_fraction) < tol])
    return r + 1j * im
# फॉरवर्ड/इनवर्स डिस्क्रीट फ़ूरिएर सीरीज़ (स्केलिंग में समस्या)
def fs_transform_bad(data, halfspan, inv=False):
    win_len = 2 * halfspan + 1  # used both as array length and frequency scale
    out = np.zeros(win_len, dtype=complex)
    if inv:
        for n in range(-halfspan, halfspan + 1):
            for k in range(-halfspan, halfspan + 1):
                out[n + halfspan] += data[k + halfspan] * cmath.exp(JJ * 2 * PI * k * n / win_len)
    else:
        for k in range(-halfspan, halfspan + 1):
            for n in range(-halfspan, halfspan + 1):
                out[k + halfspan] += (1 / win_len) * data[n + halfspan] * cmath.exp(-JJ * 2 * PI * k * n / win_len)
    return out
# Spectrum of rect
spec_a = fs_transform_bad(rect_sig, L)
spec_a = snap_near_integers(spec_a)
# Stretch in frequency by inserting zeros and scaling
factor = 5
spec_f = np.zeros(factor * M, dtype=complex)
spec_f[::factor] = 0.2 * spec_a
# Reconstruct with the same transform
guess_L = (len(spec_f) - 1) // 2
time_f = fs_transform_bad(spec_f, guess_L, inv=True)
time_f = snap_near_integers(time_f)

गलती कहाँ हुई और क्यों होती है

मूल कारण D और N के बीच इंडेक्सिंग का असंगति है, जिससे DC बिन की स्थिति बिगड़ जाती है और स्पेक्ट्रम खिसक जाता है। यदि मूल आधा‑रेंज D=L है और लंबाई N=2D+1 है, तो N को सीधे किसी फैक्टर से गुणा करने पर नई लंबाई कभी सम भी हो सकती है या, विषम होने पर भी, “बीच” का इंडेक्स DC से हट सकता है। इस उदाहरण में factor 5 के साथ ऐरे में मध्य तो है, लेकिन स्ट्राइड से ज़ीरो डालते समय DC टर्म ठीक केंद्र पर नहीं आता। जो शून्य‑आवृत्ति वाला गुणांक पहले इंडेक्स L पर था, वह कुछ इंडेक्स खिसक जाता है; यही शिफ्ट रिकंस्ट्रक्शन में आता है और काल्पनिक घटक पैदा करता है।

इसे सीधे पैड किए गए स्पेक्ट्रम को जाँचकर देखा जा सकता है। खिंचे हुए ऐरे का केंद्र इंडेक्स वहाँ नहीं है जहाँ मूल का DC सैंपल ज़ीरो डालने के बाद पहुँचता है, इसलिए इनवर्स ट्रांसफॉर्म एक ग़लत‑सँरेखित DC बिन का उपयोग करता है। सिर्फ़ यह मिसअलाइनमेंट ही फेज़ ऑफ़सेट और अप्रत्याशित इमैजिनरी भाग को समझा देता है, जबकि मूल a[n] पूरी तरह वास्तविक और सम था।

एक दूसरा, क़रीबी मुद्दा खुद डिस्क्रीट फ़ूरिएर सीरीज़ के इम्प्लीमेंटेशन में है। ऐरे की लंबाई 2D+1 है, लेकिन एक्सपोनेंट के भीतर उपयोग किया जाने वाला फ़्रीक्वेंसी स्केल 2D पर आधारित होना चाहिए, 2D+1 पर नहीं। 2D+1 को एक साथ ऐरे की लंबाई और कोण के हर में लगाना प्रभावी ग्रिड को बदल देता है और बाद में स्ट्रेच करने पर शिफ्ट और बढ़ जाता है। ऐरे का आकार भले 2D+1 रहे, लेकिन एक्सपोनेंट में नॉर्मलाइज़ेशन फ़ैक्टर 2D होना चाहिए।

उपाय: N नहीं, D को फैलाएँ, और सही फ़्रीक्वेंसी स्केल लें

दो सटीक बदलाव समस्या सुलझा देते हैं और सिग्नल की पाँच वास्तविक, ठीक‑ठाक संरेखित प्रतियाँ मिलती हैं। पहला, N = 2D + 1 को अपरिवर्तनीय मान कर रखें। स्ट्रेच करते समय D को फैक्टर से गुणा करें और फिर उसी D से N पुनर्निर्मित करें। दूसरा, फ़ूरिएर सीरीज़ वाले कोड में, जब तक ऐरे में 2D+1 सैंपल हैं, एक्सपोनेंट के नॉर्मलाइज़ेशन के लिए 2D का उपयोग करें। इन सुधारों के बाद DC गुणांक केंद्र इंडेक्स पर आता है और इनवर्स ट्रांसफॉर्म अपेक्षित वास्तविक पैटर्न लौटाता है।

import numpy as np
import cmath
L = 1000
JJ = complex(0, 1)
PI = np.pi
M = 2 * L + 1
# Rect सिग्नल
rect_sig = np.zeros(M)
for t in range(-99, 100):
    rect_sig[t + L] = 1
# FP सफ़ाई
eps_cut = 1e-10
def snap_near_integers(arr, tol=eps_cut):
    r = np.real(arr)
    im = np.imag(arr)
    r[np.abs(r) < tol] = 0
    im[np.abs(im) < tol] = 0
    r_fraction = r - np.round(r)
    im_fraction = im - np.round(im)
    r[np.abs(r_fraction) < tol] = np.round(r[np.abs(r_fraction) < tol])
    im[np.abs(im_fraction) < tol] = np.round(im[np.abs(im_fraction) < tol])
    return r + 1j * im
# सही डिस्क्रीट फ़ूरिएर सीरीज़: ऐरे आकार 2D+1, फ़्रीक्वेंसी स्केल 2D
def fs_transform_ok(data, halfspan, inv=False):
    freq_span = 2 * halfspan            # scale for the exponent
    out = np.zeros(freq_span + 1, dtype=complex)  # array length still 2D+1
    if inv:
        for n in range(-halfspan, halfspan + 1):
            for k in range(-halfspan, halfspan + 1):
                out[n + halfspan] += data[k + halfspan] * cmath.exp(JJ * 2 * PI * k * n / freq_span)
    else:
        for k in range(-halfspan, halfspan + 1):
            for n in range(-halfspan, halfspan + 1):
                out[k + halfspan] += (1 / freq_span) * data[n + halfspan] * cmath.exp(-JJ * 2 * PI * k * n / freq_span)
    return out
# फॉरवर्ड ट्रांसफॉर्म
spec_a = fs_transform_ok(rect_sig, L)
spec_a = snap_near_integers(spec_a)
# गुणांकों के बीच शून्य डालकर स्ट्रेच करना, फैक्टर से स्केल करना
factor = 5
L_new = factor * L
M_new = 2 * L_new + 1
spec_g = np.zeros(M_new, dtype=complex)
spec_g[::factor] = spec_a / factor
# इनवर्स ट्रांसफॉर्म
time_g = fs_transform_ok(spec_g, L_new, inv=True)
time_g = snap_near_integers(time_g)

यदि आप टूटे हुए संस्करण में DC‑बिन का मिसअलाइनमेंट देखना चाहें, तो एक छोटा‑सा प्रॉब इसे साफ़ कर देता है। केंद्र एंट्री ज़ीरो इन्सर्शन के बाद DC मान नहीं रहती; इस उदाहरण (factor 5) में वह दो इंडेक्स इधर‑उधर बैठी मिलती है।

# ग़लत तरीक़े में DC शिफ्ट दिखाना
factor = 5
spec_f = np.zeros(factor * M, dtype=complex)
spec_f[::factor] = 0.2 * spec_a
mid = spec_f.shape[0] // 2
print(spec_f[mid])       # ends up 0
print(spec_f[mid - 2])   # holds the scaled DC term

यह क्यों महत्वपूर्ण है

जब आप फ़ूरिएर सीरीज़ को इंडेक्स के आधार पर बदलते हैं, तो आधा‑रेंज D, लंबाई N = 2D + 1 और DC बिन का संबंध महज़ बहीखाता नहीं है; यही वह ग्रिड तय करता है जिस पर फ़्रीक्वेंसी ऑपरेशन होते हैं। सीधे N को छेड़कर स्ट्रेच करना और एक्सपोनेंशियल के स्केल को गड्ड‑मड्ड करना पूरे स्पेक्ट्रम को ऑफ़सेट कर देता है, इसलिए इनवर्स ट्रांसफॉर्म को एक खिसके हुए फ़्रीक्वेंसी ग्रिड से पुनर्निर्माण करना पड़ता है। यही अप्रत्याशित इमैजिनरी घटक और समय क्षेत्र में साफ़‑सुथरी आवधिक प्रतियाँ खोने की वजह है।

मुख्य बातें

N = 2D + 1 को सख़्त नियम की तरह बनाए रखें और N नहीं, D को फैलाएँ। डिस्क्रीट फ़ूरिएर सीरीज़ में, जब सिग्नल में अभी भी 2D + 1 सैंपल हैं, एक्सपोनेंट के नॉर्मलाइज़ेशन के लिए 2D का उपयोग करें। स्पेक्ट्रल सैंपलों के बीच शून्य जोड़ते समय यह सुनिश्चित करें कि DC टर्म नए स्पेक्ट्रम के सही केंद्र इंडेक्स पर पहुँचे। ये तीन काम कर दें, तो इनवर्स अपेक्षित वास्तविक, दोहराई गई वेवफ़ॉर्म देगा, बिना अनचाहे फेज़ आर्टिफ़ैक्ट के। अगर प्रदर्शन चिंता बनता है, तो भीतर के समेशन को वेक्टराइज़ करने पर विचार करें, लेकिन यहाँ शुद्धता पूरी तरह इंडेक्सिंग और फ़्रीक्वेंसी स्केल की है।

यह लेख StackOverflow पर प्रश्न (लेखक: Nate3384) और chrslg के उत्तर पर आधारित है।