1

私は現在、次の本を研究しています：VidiSaptariによる「FourierTransformSpectroscopyInstrumentationEngineering」。私の質問は、本の付録Cのコードに基づいて、以下のコードに関連しています。以下のコードは、波数[cm-1] 5000、10000、および15000の3つの波のインターフェログラムを計算し、FFTを実行します。情報を取得します。スケーリングされていない出力の大きさは、1ではなく1600です。

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になります。

同じリンクの2番目の解に従い、 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 (したがって電力 0.5) の単一 cos 波を表しています。

個人的には、MATLAB の FFT 出力を でスケーリングし1/sqrt(N)、その IFFT 出力を でスケーリングすることを好みますsqrt(N)。このようにして、時間領域と周波数領域のシーケンスのエネルギーは常に等しくなります。

于 2012-08-31T07:55:54.303 に答える
0

if you want the same energy in the IFFT input and ouput you have to multiply the IFFT output (time signal) by sqrt(N) where N is the size of the transform.

here's the code :

same thing for the FFT (divide output by sqrt(N)); hope it helps 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

于 2013-07-31T21:49:52.383 に答える