FFT归一化

Question

我知道这个问题被问到令人作呕，但不知何故我无法让它正常工作。我创建了一个具有单位幅度的 440 Hz 正弦波。现在，在 FFT 之后，440 Hz 的 bin 有一个明显的峰值，但该值不正确。我希望看到 0 dB，因为我正在处理单位幅度的正弦波。相反，计算的功率远高于 0 dB。我使用的公式很简单

for (int i = 0; i < N/2; i++) 
{  
    mag = sqrt((Real[i]*Real[i] + Img[i]*Img[i])/(N*0.54)); //0.54 correction for a Hamming Window

    Mag[i] = 10 * log(mag) ;      
}

我可能应该指出，时域中的总能量等于频域中的能量（Parseval 定理），所以我知道我的 FFT 例程很好。

非常感谢任何帮助。

Answer 1

我一直在为工作而苦苦挣扎。看来很多软件例程/书籍对FFT的归一化有点草率。 我得到的最好的总结是：能量需要守恒——这是 Parseval 的定理。此外，在 Python 中编码时，您很容易丢失一个元素并且不知道它。请注意，numpy.arrays 索引不包括最后一个元素。

a = [1,2,3,4,5,6]
a[1:-1] = [2,3,4,5]
a[-1] = 6

这是我正确规范化 FFT 的代码：

# FFT normalization to conserve power
    
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal
    
    
sample_rate = 500.0e6
time_step = 1/sample_rate
    
carrier_freq = 100.0e6
    
    
# number of digital samples to simulate
num_samples = 2**18 # 262144
t_stop = num_samples*time_step
t = np.arange(0, t_stop, time_step)
    
# let the signal be a voltage waveform,
# so there is no zero padding
carrier_I = np.sin(2*np.pi*carrier_freq*t)
    
#######################################################
# FFT using Welch method
# windows = np.ones(nfft) - no windowing
# if windows = 'hamming', etc.. this function will
# normalize to an equivalent noise bandwidth (ENBW)
#######################################################
nfft = num_samples  # fft size same as signal size 
f,Pxx_den = scipy.signal.welch(carrier_I, fs = 1/time_step,\
                        window = np.ones(nfft),\
                        nperseg = nfft,\
                        scaling='density')
    
#######################################################
# FFT comparison    
#######################################################
    
integration_time = nfft*time_step
power_time_domain = sum((np.abs(carrier_I)**2)*time_step)/integration_time
print 'power time domain = %f' % power_time_domain
    
# Take FFT.  Note that the factor of 1/nfft is sometimes omitted in some 
# references and software packages.
# By proving Parseval's theorem (conservation of energy) we can find out the 
# proper normalization.
signal = carrier_I 
xdft = scipy.fftpack.fft(signal, nfft)/nfft 
# fft coefficients need to be scaled by fft size
# equivalent to scaling over frequency bins
# total power in frequency domain should equal total power in time domain
power_freq_domain = sum(np.abs(xdft)**2)
print 'power frequency domain = %f' % power_freq_domain
# Energy is conserved
    
    
# FFT symmetry
plt.figure(0)
plt.subplot(2,1,1)
plt.plot(np.abs(xdft)) # symmetric in amplitude
plt.title('magnitude')
plt.subplot(2,1,2) 
plt.plot(np.angle(xdft)) # pi phase shift out of phase
plt.title('phase')
plt.show()
    
xdft_short = xdft[0:nfft/2+1] # take only positive frequency terms, other half identical
# xdft[0] is the dc term
# xdft[nfft/2] is the Nyquist term, note that Python 2.X indexing does NOT 
# include the last element, therefore we need to use 0:nfft/2+1 to have an array
# that is from 0 to nfft/2
# xdft[nfft/2-x] = conjugate(xdft[nfft/2+x])
Pxx = (np.abs(xdft_short))**2 # power ~ voltage squared, power in each bin.
Pxx_density = Pxx / (sample_rate/nfft)  # power is energy over -fs/2 to fs/2, with nfft bins
Pxx_density[1:-1] = 2*Pxx_density[1:-1] # conserve power since we threw away 1/2 the spectrum
# note that DC (0 frequency) and Nyquist term only appear once, we don't double those.
# Note that Python 2.X array indexing is not inclusive of the last element.
    
    
    
Pxx_density_dB = 10*np.log10(Pxx_density)
freq = np.linspace(0,sample_rate/2,nfft/2+1) 
# frequency range of the fft spans from DC (0 Hz) to 
# Nyquist (Fs/2).
# the resolution of the FFT is 1/t_stop
# dft of size nfft will give nfft points at frequencies
# (1/stop) to (nfft/2)*(1/t_stop)
    
plt.figure(1)
plt.plot(freq,Pxx_density_dB,'^')
plt.figure(1)
plt.plot(f, 10.0*np.log10(Pxx_den))
plt.xlabel('Freq (Hz)'),plt.ylabel('dBm/Hz'),#plt.show()
plt.ylim([-200, 0])
plt.show()

Answer 2

许多常见（但不是全部）FFT 库将单位幅度正弦曲线的 FFT 结果按 FFT 的长度缩放。这保持了 Parsevals 等式，因为较长的正弦曲线比相同幅度的较短的正弦曲线代表更多的总能量。

如果您不想在使用这些库之一时按 FFT 长度进行缩放，则在计算以 dB 为单位的幅度之前除以长度。

Answer 3

归一化可以通过许多不同的方式完成 - 取决于窗口、样本数量等。

常用技巧：对已知信号进行 FFT 并按峰值进行归一化。假设在上面的示例中，您的峰值是 123 - 如果您希望它是 1，则将其（以及使用此算法获得的所有结果）除以 123。