需要帮助:FFT输出缩放

Question

我目前正在学习Vidi Saptari的《傅里叶变换光谱仪器工程》一书。我的问题与以下代码有关，该代码基于书附录C中的代码。以下代码计算波数分别为[cm-1]5000、10000和15000的3个波的干涉图，然后执行FFT以检索信息。未缩放的输出的幅值为1600，而不是1。

clear;
% Sampling clock signal generation
samp_period_nm = 632.8 / 4; % sampling period in nm. 632.8 is the HeNe laser's wavelength
samp_period = 1 * samp_period_nm * 10^-7; % sampling period in cm. 
scan_dist   = 0.1; % mirror scan distance in cm. 
no_elements = floor(scan_dist/samp_period);
x_samp = 0:samp_period:samp_period*no_elements; %Vector of clock signals in cm
xn_samp = x_samp .* (1 + rand(1, length(x_samp)));
v1 = 5000;
v2 = 10000;
v3 = 15000;
arg  = 4 * pi * x_samp;
y   = cos(arg*v1) + cos(arg*v2) + cos(arg*v3) ;
total_data = 2^18;
no_zero_fills=[total_data - length(y)];
zero_fills=zeros(1, no_zero_fills);
%triangular apodization
n_y = length(y);
points = 1:1:n_y;
tri = 1 - 1/(n_y) * points(1:n_y);
y = y.*tri; %dot product of interferogram with triangular apodization function
y = [y zero_fills];   %zero filling
% FFT operation
fft_y = fft(y);
% fft_y = fft_y / n_y;
% fft_y = fft_y * samp_period;
fft_y(1) = [];
n_fft=length(fft_y);
spec_y = abs(fft_y(1:n_fft/2)); %spectrum generation
nyquist = 1 / (samp_period * 4);
freq = (1:n_fft/2)/(n_fft/2)*nyquist; %frequency scale generation
figure();
plot(freq, spec_y);   % plot of spectrum vs wave number
xlabel('Wavenumber [cm-1]');
ylabel('Intesity [V]');

通过将fft（fft_y）的结果与dt＝samp_period相乘，如这里所建议的，峰值为0.025。

根据同一链路的第二个解，通过将fft_y除以n_y（y的长度），大小为0.25。

很明显，我做错了什么。感谢您的帮助。

谢谢，

Answer 1

这里唯一做错的是期望光谱中的峰值为1。根据DFT的Parseval定理，时域信号的能量等于频域信号的能量除以序列N的长度

td_energy = sum( abs(y).^2 )
fd_energy = sum( abs(fft_y).^2 )
td_energy - fd_energy / length(y) % won't be exactly zero because you deleted the zero frequency bin.

因此，频谱中的峰值并不代表余弦波在时域中的振幅，而是它们的能量。还要注意，在这一点上，能量比你可能预期的要低，因为你填充了很多零。

在实践中，特定频率的平均功率通常更令人感兴趣。考虑以下代码示例

t = linspace(-4*pi, 4*pi, 2^16);
N = length(t); % DFT length
y = cos(t); % single cosine wave
y_pow = sum( abs(y).^2 ) / N; % is 0.5
fft_y = fft(y);
fft_y_pow = (sum( abs(fft_y).^2 ) / N) /N; % is 0.5
figure; plot(abs(fft_y)./N);

功率是通过将能量除以序列N的长度来获得的。如果将频谱除以N，则获得每个频率的平均功率。在上面的例子中，你可以识别出一个高度为0.5的单个峰值，它代表振幅为1（因此功率为0.5）的单个cos波

就我个人而言，我更喜欢用1/sqrt(N)缩放MATLAB的FFT输出，用sqrt(N)缩放其IFFT输出。这样，时域序列和频域序列的能量总是相等的。

Answer 2

如果你想在IFFT输入和输出中获得相同的能量，你必须相乘IFFT输出（时间信号）乘以sqrt（N），其中N是变换的大小。

这是代码：

FFT也是如此（将输出除以sqrt（N））；希望能有所帮助Felix

N = 4096;
Freq = randn(N,1)+1j*randn(N,1);
Time = sqrt(N)*ifft(Freq,N);
FreqEn = sum(real(Freq).^2 + imag(Freq).^2); 
TimeEn = sum(real(Time).^2 + imag(Time).^2); 

TimeEn/Freq