1

我目前正在研究以下书籍：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。

显然，我做错了什么。任何帮助表示赞赏。

谢谢，

4

2 回答 2

2

您在这里做错的唯一一件事是期望频谱中的峰值为 1。根据 Parseval 的 DFT 定理，时域信号的能量等于频域信号的能量除以序列的长度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 的单个 cos 波（因此功率为 0.5）。

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

于 2012-08-31T07:55:54.303 回答
0

如果您希望 IFFT 输入和输出中的能量相同，则必须将 IFFT 输出（时间信号）乘以 sqrt(N)，其中 N 是变换的大小。

这是代码：

FFT 也是一样的（输出除以 sqrt(N)）；希望它可以帮助菲利克斯

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);

时间/频率

于 2013-07-31T21:49:52.383 回答