信息处理 - 不带数字滤波器的加权信号 - 吾爱随笔录

不带数字滤波器的加权信号

信息处理 fft Python

2022-01-29 17:47:20

我本质上想通过fft在 Python 中使用对信号进行加权。但我想我犯了一些逻辑错误。

我的步骤：

初始化正弦波：y = np.array(np.rint(32768 * np.sin(1000 * 2.0 * np.pi * np.linspace(0,1,44100))))
取 RMS：print(20.0 *np.log10(rms_flat(y)))- 给我 $87\textrm{ dB}$

fft带汉宁窗：

yf = np.abs(np.fft.fft(samples * np.hanning(len(samples))))

取 RMS：给我大约 - 为什么不再是？ $130\textrm{ dB}$ $87\textrm{ dB}$
最后是权重谱 - 见下文 - 获得 Inf RMS ......

将 matplotlib.pyplot 导入为 plt
将 numpy 导入为 np


#samplerate
Fs = 44100

#samplecount
A = Fs

#Sample Interval
Ts = 1/Fs

#Zeitvektor
t = np.linspace(0.0, Fs*Ts, Fs)

#Sinus Freq
ff = 1000

y = np.array(np.rint(32768 * np.sin(ff * 2.0 * np.pi * t)))

fig = plt.figure(figsize=(10,10))
fig.subplots_adjust(hspace=0.7)
fig.subplots_adjust(wspace=0.25)
plt.subplot(321)
plt.title("Samples")
plt.ylabel("Amplitude")
plt.xlabel("Time")
plt.plot(t, y)


def rms_flat(a):  # from matplotlib.mlab
    """
    Return the root mean square of all the elements of *a*, flattened out.
    """
    return np.sqrt(np.mean(np.absolute(a)**2))


print(20.0 *np.log10(rms_flat(y)))


def fft(samples):
    yf = np.abs(np.fft.fft(samples * np.hanning(len(samples))))
    xf = np.linspace(0.0, 1 / (2.0 * Ts), Fs / 2)
    plt.subplot(322)
    plt.title("Abs Pos Spectrum")
    plt.ylabel("Amplitude")
    plt.xlabel("Frequency")
    plt.plot(xf, 2.0/A * yf[:A//2])
    print(20.0 * np.log10(rms_flat(yf)))
    return(yf, xf)

yf, xf = fft(y)


def aWeighting(yf, xf):
    f1 = 20
    f2 = 107
    f3 = 737
    f4 = 12194
    N = len(xf)
    i = 0
    yfa = []
    while i < N:
        a = 20.0 * np.log10((f4 ** 2 * xf[i] ** 4) / ((xf[i] ** 2 + f1 ** 2) * np.sqrt(xf[i] ** 2 + f2 ** 2) * np.sqrt(xf[i] ** 2 + f3 ** 2) * (xf[i] ** 2 + f4 ** 2))) + 2.0
        yfa.append(20.0 * np.log10(yf[i]) + a)
        i += 1

    yfa = np.array(yfa)
    print(rms_flat(yfa))
    plt.subplot(323)
    plt.title("A-weighted Spektrum in dB")
    plt.ylabel("dB")
    plt.xlabel("Frequency")
    plt.plot(xf, 2.0/A * yfa[:A//2])
    plt.show()
    return(yfa, weight)

yfa, weight = aWeighting(yf, xf)

1个回答

FFT 的处理增益为 10Log(N)，适用于标准偏差的测量，其中 N 是信号的长度，汉明窗的相干增益为 0.54（这是一个损失）。请参阅 Fred harris 的著名论文On the Use of Windowing，其中详细介绍了这些收益和损失。

让我们看看这是否可行。您有 44100 个样本，因此净处理增益为 $10Log10(44100\times.54) = 43.8 \text{ dB}$

您的峰值正弦波为 32768，因此 rms 应为，这与您的测量结果一致。 $20log_{10}(32768/\sqrt{2})=87.3 \text { dB}$

获取窗口和 FFT 后，您看到 RMS 信号增加到 130 dB，而在 FFT 之前测量到 87 dB。130-87 = 43 分贝。

43 dB 的差异与上面给出的净处理增益公式所预测的一致。

以下是更多详细信息，可能有助于进一步了解这是如何发生的：

输入信号的相干增益为 N。输入信号的时间为

A s i n (ω t) = \frac{A}{2} \frac{e^{j ω t}}{j} - \frac{A}{2} \frac{e^{j ω t}}{j}

$A sin(\omega t) = \frac{A}{2}\frac{e^{j\omega t}}{j}-\frac{A}{2}\frac{e^{j\omega t}}{j}$

从上面的欧拉展开式中看到的重要项目是正弦由两个指数组成，每个指数的幅度为正弦波本身的一半。在 FFT 之后，由于 FFT 中的相干处理增益，这些中的每一个都增加了 N 倍。

因此，FFT 中的两个音调中的每一个都具有的幅度（如果您以精确的整数个样本进行连贯采样，这将是准确的，否则我们会看到它略有减少，因为“扇形损失”）。 $\frac{AN}{2}$

如果连贯地采样，除了所描述的两个音调之外，其他每个 bin 都将为零（我们可以通过 Parseval 定理进一步展示信号中所有能量守恒的一般情况，但为了简单起见，我将仅描述相干的情况采样，意味着采样时钟是正弦波的整数倍）。

因此在计算 rms 时，由于只有 2 个非零值但我们使用 N 个样本，FFT 的 rms 计算将减少为：

σ = \sqrt{\frac{{(\frac{A N}{2})}^{2} + {(\frac{A N}{2})}^{2}}{N}}

$\sigma = \sqrt{\frac{\left(\frac{AN}{2}\right)^2+\left(\frac{AN}{2}\right)^2}{N}}$

这减少到

σ = \sqrt{N} \frac{A}{\sqrt{2}}

$\sigma = \sqrt{N}\frac{A}{\sqrt{2}}$

我们将其识别为乘以输入正弦波的 rms。当我们在进行 FFT 之前对信号进行窗口化时，相干幅度（先前增加 N）现在将被窗口的相干增益降低，对于汉明窗口的情况，该增益为 0.54。 $\sqrt{N}$

因此，对于使用汉明窗的情况，FFT 的标准偏差预计倍。所以以分贝为单位： $\sqrt{(N\times0.54)}$

20 l o g_{10} (σ) = 10 l o g_{10} (N \times 0.54)

$20log_{10}(\sigma) = 10log_{10}(N\times0.54)$

其它你可能感兴趣的问题

上一篇在一维 DCT 中，为什么输入是向量？下一篇如何在 Matlab 中实现 Viola-Jones 算法？