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结果。由于较长的正弦曲线比相同振幅的较短正弦曲线代表更多的总能量，因此这保持了部分相等。

如果您在使用这些库时不想按FFT长度缩放，那么在计算dB中的幅值之前，请除以长度。

Answer 3

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

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