IAAFT के साथ log‑normal सिग्नल: blue/violet शोर के लक्ष्य PSD को बनाए रखने की गाइड

Exponentiation से blue/violet शोर का PSD क्यों ढहता है? IAAFT से log‑normal marginal और target PSD साथ‑साथ मिलाएँ: log–log इंटरपोलेशन, कोड और सलाह.

बढ़ती power spectral density वाले लॉग‑नॉर्मल सिग्नल बनाना कागज़ पर आसान लगता है: Gaussian शोर लें, उसका स्पेक्ट्रम आकार दें, और उसे exponentiate कर दें। व्यवहार में, आखिरी कदम उच्च‑पास प्रकारों—जैसे blue या violet शोर—की स्पेक्ट्रल संरचना बिगाड़ देता है। यह मार्गदर्शिका बताती है कि ऐसा क्यों होता है और IAAFT की मदद से लक्ष्य marginal वितरण (उदा., log‑normal) और लक्ष्य PSD—दोनों से मेल खाते सिग्नल कैसे तैयार करें।

समस्या का संक्षेप

उद्देश्य ऐसे रैंडम अनुक्रम बनाना है जिनके नमूने धनात्मक मानों पर समर्थित, लंबी पूंछ वाले वितरण का पालन करें (उदा., log‑normal), और जिनकी power spectral density कम आवृत्तियों पर कम हो तथा आवृत्ति बढ़ने पर बढ़े (उदा., blue शोर +10 dB/decade या violet शोर +20 dB/decade)। Log‑normal वितरण से सैंपल लेना आसान है। FFT के जरिए Gaussian शोर को लक्ष्य PSD के अनुरूप shape करना भी सामान्य अभ्यास है। मुश्किल तब आती है जब पहले से shape किए गए प्रोसेस पर एक nonlinear mapping (exponential) लागू कर दी जाती है।

पुनरुत्पाद्य उदाहरण: पहले shaping, फिर exponentiation

नीचे दिया कोड कई मानक PSDs के साथ colored Gaussian शोर बनाता है, फिर log‑normal नमूने पाने के लिए exponential mapping लगाता है। Exponentiation के बाद white और pink शोर के लिए PSD का आकार बना रहता है, लेकिन blue और violet शोर का स्पेक्ट्रम सफेद‑जैसा हो जाता है।

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch
# rFFT के जरिए spectral shaping: target magnitude स्पेक्ट्रम से गुणा करें
# और inverse लें.
def shape_spectrum(n_points, psd_fn=lambda f: 1.0):
    freq_side = np.fft.rfft(np.random.randn(n_points))
    spec = psd_fn(np.fft.rfftfreq(n_points))
    spec = spec / np.sqrt(np.mean(spec ** 2))
    colored = freq_side * spec
    return np.fft.irfft(colored)
# सहायक: PSD फ़ंक्शन को length->signal जनरेटर में बदलने के लिए.
def with_psd(psd_fn):
    return lambda n: shape_spectrum(n, psd_fn)
@with_psd
def make_white(f):
    return 1.0
@with_psd
def make_blue(f):
    return np.sqrt(f)
@with_psd
def make_violet(f):
    return f
@with_psd
def make_brown(f):
    return 1.0 / np.where(f == 0, float('inf'), f)
@with_psd
def make_pink(f):
    return 1.0 / np.where(f == 0, float('inf'), np.sqrt(f))
def synthesize_colored(n_points, kind="white"):
    if kind == "white":
        return make_white(n_points)
    elif kind == "pink":
        return make_pink(n_points)
    elif kind == "blue":
        return make_blue(n_points)
    elif kind == "violet":
        return make_violet(n_points)
    elif kind == "brown":
        return make_brown(n_points)
# डेमो: विभिन्न PSD प्रकारों पर Gaussian बनाम log‑normalized संस्करणों की तुलना.
n_points = 100_000
kinds = ["white", "pink", "blue", "violet"]
colors = ["lightgray", "magenta", "blue", "violet"]
fig, axes = plt.subplots(2, 2, sharex=False, sharey="row")
bins_gauss = np.linspace(-5, 5, 101)
for nm, col in zip(kinds, colors):
    x = synthesize_colored(n_points, nm)
    axes[0, 0].hist(x, bins=bins_gauss, color=col, label=nm, histtype="step", density=True)
    fx, px = welch(x)
    axes[1, 0].loglog(fx, px, color=col, label=nm)
bins_logn = np.linspace(0, 10, 101)
for nm, col in zip(kinds, colors):
    y = np.exp(synthesize_colored(n_points, nm))
    axes[0, 1].hist(y, bins=bins_logn, color=col, label=nm, histtype="step", density=True)
    fy, py = welch(y)
    axes[1, 1].loglog(fy, py, color=col, label=nm)
axes[0, 0].set_title("Gaussian samples")
axes[0, 1].set_title("Exp(Gaussian) samples")
axes[0, 0].set_ylabel("Density")
axes[0, 0].set_xlabel("Value")
axes[0, 1].set_xlabel("Value")
axes[1, 0].set_ylabel("Power")
axes[1, 0].set_xlabel("Frequency [Hz]")
axes[1, 1].set_xlabel("Frequency [Hz]")
axes[0, 0].legend(loc="upper left")
fig.tight_layout()
plt.show()

Exponentiation के बाद स्पेक्ट्रल आकार क्यों ढह जाता है

एक्सपोनेन्शिएशन एक nonlinear क्रिया है। रैखिक फ़िल्टरिंग के विपरीत, यह उस रैखिक सहसंबंध संरचना को सुरक्षित नहीं रखती जिसे PSD निरूपित करता है। High‑pass स्पेक्ट्रा (जैसे blue और violet शोर) के मामले में, exponential लगाने से स्पेक्ट्रम काफी बदल जाता है और परिणामी PSD सफेद शोर जैसा दिखने लगता है। यही देखा गया व्यवहार है।

एक अलग, स्वतंत्र बात यह है कि औसत ऊर्जा बनाए रखने के लिए shaping कर्व को normalize करना—जैसा कि फ़्रीक्वेंसी डोमेन में colored शोर बनाते समय अक्सर किया जाता है—बाद की nonlinear विकृति की भरपाई नहीं करता। Whitening प्रभाव normalization से नहीं, बल्कि nonlinear mapping से आता है।

IAAFT: PSD और marginal वितरण को एक साथ लागू करना

दोनों बाधाओं—लक्ष्य marginal वितरण और लक्ष्य PSD—को साथ में पूरा करने का एक ठोस तरीका है उनके बीच आवर्तन करना। Iterative Amplitude Adjusted Fourier Transform (IAAFT) एल्गोरिदम यही करता है। यह दो प्रोजेक्शन बारी‑बारी से लागू करता है: एक चरण में वांछित magnitude स्पेक्ट्रम लगाया जाता है; दूसरे में रैंक‑रीमैपिंग द्वारा लक्ष्य नमूना वितरण बहाल किया जाता है। यह तरीका stochastic सिग्नलों के लिए वर्णित है और संदर्भित साहित्य में अच्छी तरह प्रलेखित है।

एक व्यावहारिक विवरण जो स्पेक्ट्रल लक्ष्य की निष्ठा को काफी सुधारता है यह है कि जब PSD को log–log रूप में दिया गया हो, तो DFT magnitude को log–log स्पेस में piecewise‑linear इंटरपोलेशन से बनाया जाए। यह उन ढलानों की परिभाषा के तरीके से मेल खाता है (जैसे “+10 dB/decade”) और PSD नियंत्रण को कहीं अधिक सटीक बनाता है।

import numpy as np
from scipy import signal as sig
from scipy.stats import rankdata
from numpy.fft import rfft, irfft
# विन्यास
n_pts = 2 ** 16
fs_rate = 1.0
max_iter = 50
eps_tol = 1e-6  # सापेक्ष PSD त्रुटि सीमा
# लक्ष्य marginal: Gaussian ड्रॉज़ को exponentiate करके log‑normal आकार
target_seed = np.exp(np.random.randn(n_pts))
asc_template = np.sort(target_seed)
# लक्ष्य magnitude स्पेक्ट्रम: उदाहरण +10 dB/decade ("blue noise")
# log–log अक्षों पर piecewise‑linear इंटरपोलेशन से बनाएं.
brk_freqs = np.array([1e-5 * fs_rate, fs_rate / 2])
brk_psd_db = 10 * np.log10(brk_freqs)
freq_bins = np.fft.rfftfreq(n_pts, d=1 / fs_rate)
# log10(0) से बचने को पहला बिंदु जोड़ें और DC को स्पष्ट रूप से सेट करें.
brk_freqs = np.insert(brk_freqs, 0, fs_rate / n_pts)
brk_psd_db = np.insert(brk_psd_db, 0, brk_psd_db[0])
psd_interp_db = np.interp(
    np.log10(freq_bins[1:]),
    np.log10(brk_freqs[1:]),
    brk_psd_db[1:]
)
# amplitude स्पेक्ट्रम में बदलें; high‑pass रूपों के लिए DC बिन को शून्य रखें.
target_amp = np.empty_like(freq_bins)
target_amp[0] = 0.0
target_amp[1:] = 10 ** (psd_interp_db / 20)
# रैंक‑आधारित रीमैपिंग: समान रैंक वाले लक्ष्य marginal के मानों से प्रतिस्थापित करें.
def remap_by_rank(values, sorted_ref):
    ranks = rankdata(values, method='ordinal') - 1
    return sorted_ref[ranks]
# IAAFT लूप: स्पेक्ट्रम लागू करना और marginal पुनर्निर्माण बारी‑बारी से.
synth = target_seed.copy()
prev_err = np.inf
for i in range(max_iter):
    # लक्ष्य magnitude स्पेक्ट्रम लागू करें.
    X = rfft(synth)
    X = target_amp * np.exp(1j * np.angle(X))
    shaped = irfft(X, n=n_pts)
    # रैंक‑रीमैपिंग से लक्ष्य marginal बहाल करें.
    synth = remap_by_rank(shaped, asc_template)
    # समय‑समय पर अभिसरण जाँचें.
    if (i % 5 == 0) or (i == max_iter - 1):
        cur_err = np.linalg.norm(np.abs(rfft(synth)) - target_amp) / np.linalg.norm(target_amp)
        print(f"iter {i:2d}: PSD rel-err = {cur_err:.3e}")
        if cur_err > prev_err - eps_tol:
            print("converged.")
            break
        prev_err = cur_err

फ़िल्टरिंग‑आधारित तरीकों की तुलना

Blue‑noise shaping फ़िल्टर्स (FIR या IIR) के साथ log‑normal नमूनों की रैखिक फ़िल्टरिंग वांछित PSD ढलानें दे सकती है। परिणाम सख्ती से धनात्मक और दाएँ तरफ झुका हुआ रहता है, पर नमूनों के रैखिक मिश्रण के कारण वह log‑normal नहीं रहता। उपलब्ध परिणामों में, IAAFT ने परखे गए तरीकों में PSD और marginal वितरण—दोनों का एक साथ सबसे अच्छा मेल दिया, और साथ ही दशकों में परिभाषित PSD लक्ष्यों को ट्रैक करने के लिए log–log इंटरपोलेशन का लाभ मिला।

फिर भी, परीक्षणों में पाया गया कि log‑normal लक्ष्य वितरण और non‑white PSD के साथ, प्रेक्षित marginal भटक सकता है और एक की जगह शून्य पर मोड दिखा सकता है—हालाँकि पूंछ लंबी बनी रहती है। परिणामों का सत्यापन करते समय यह एक महत्वपूर्ण सावधानी है।

यह क्यों मायने रखता है

कई डाउनस्ट्रीम एल्गोरिदम विशेष PSD गुणों या marginal विशेषताओं को मान कर चलते हैं। अगर इनमें से कोई भी पक्ष बिगड़ जाए—जैसे nonlinear mapping के बाद PSD बैठ जाए, या फ़िल्टरिंग के बाद वितरण खिसक जाए—तो बेंचमार्क, डिटेक्टर या सिमुलेशन‑आधारित विश्लेषण भ्रामक हो सकते हैं। दोनों बाधाएँ साथ‑साथ लागू करने से ऐसे मिसमैच से बचाव होता है।

व्यावहारिक निष्कर्ष

यदि आपको निर्धारित PSD और कोई विशिष्ट नमूना वितरण—जैसे log‑normal—दोनों चाहिए, तो IAAFT जैसा आवर्तक ढांचा प्रभावी है। यह स्पेक्ट्रम लागू करने और रैंक‑आधारित amplitude पुनर्स्थापन के बीच अदला‑बदली करता है, और PSD निर्माण के लिए log–log इंटरपोलेशन के साथ +10 dB/decade जैसी ढलानों को बारीकी से ट्रैक कर सकता है। सत्यापन करते समय PSD और marginal वितरण—दोनों को परखें। ध्यान रखें कि non‑white लक्ष्यों में प्रेक्षित वितरण फिर भी भटक सकता है; मोड और पूंछ का व्यवहार अपेक्षाओं के विरुद्ध तो नहीं, यह जाँचें। जहाँ केवल PSD मायने रखता है वहाँ सीधी फ़िल्टरिंग सरल हो सकती है; जहाँ दोनों बाधाएँ जरूरी हों, वहाँ आवर्तक तरीका अधिक उपयुक्त है।

संक्षेप में, स्पेक्ट्रा को shape करें और amplitudes को मिलकर बहाल करें, डेटा के दोनों दृष्टिकोणों की जाँच करें, और जब PSD लक्ष्य दशकों में दिए हों तो log–log हैंडलिंग को प्राथमिकता दें।