FFT의 window energy correction 예제

목차

 윈도우는 Amplitude Sum의 크기에는 영향이 거의 없으나, Power Spectrum에는 원신호와 다른 결과를 나타낸다. 아래 예제를 통해 그 차이를 확인해 보고, Energy Correction 값을 구하는 방법도 확인해 보길 바란다.

OCTAVE 예제

% 샘플링 주파수와 신호 길이 설정
clear;
clc;
fs = 1000;  % 샘플링 주파수
N = 4096;    % 신호 길이 (FFT 포인트 수)
t = (0:N-1)/fs;

% 이상적인 사인파 신호 생성
f0 = 50; % 신호 주파수
signal = sin(2 * pi * f0 * t).*10;

% Blackman 윈도우 생성 및 적용
window = blackman(N)';
windowed_signal = signal .* window;

% FFT 계산 및 정규화
FFT_signal = fft(signal);
FFT_windowed_signal = fft(windowed_signal);
NFFT = length(FFT_signal);

% 주파수 BIN 넓이
f_bin_width = fs / NFFT;

% Power Spectrum 계산
PowerSpectrum = (abs(FFT_windowed_signal./N.*2).^2);

% 주파수 BIN 설정 (f0 근처 주파수 성분 찾기)
f_bin = round(f0 / f_bin_width);

% FFT 값의 합산 (절대값)
fft_sum = sum(abs(FFT_windowed_signal(f_bin-5:f_bin+5))) .*2/ NFFT;

% 에너지 보정 계수 계산
window_energy_correction = sqrt(sum((abs(fft(window))./N).^2) );
window_energy_correction2 = sum(abs(fft(window))./N);

% 주파수 BIN 합산 (넓게 설정, 여기서는 ±5 BIN)


power_sum = sum(PowerSpectrum(f_bin-5:f_bin+5));
power_sqrt = sqrt(power_sum);

% 보정된 Power Spectrum의 제곱근 계산
corrected_power_sqrt = power_sqrt ./ window_energy_correction;
%corrected_power_sqrt = power_sqrt;


% 결과 출력
printf('Frequency BIN width: %.3f Hz\n', f_bin_width);
printf('Square root of the sum of Power Spectrum: %.3f\n', power_sqrt);

  printf('Corrected square root of the sum of Power Spectrum: %.3f\n',
  corrected_power_sqrt);

printf('Sum of the absolute FFT values: %.3f\n', fft_sum);

윈도우 weight를 계산하는 function은 아래와 같이 생성해서 사용하면 편리하다. "fft_weight.m" 파일 예제

function weight_result = fft_weight(window, NFFT) % Calculate window weight

  window_fft=abs(fft(window,NFFT)./NFFT);
  weight_result=1/sqrt(sum(window_fft.^2));

  return;
endfunction

Test Bench "tb_fft_weight.m" 파일 예제

% 샘플링 주파수와 신호 길이 설정
clear;
clc;
addpath="E:\Document\project\OCTAVE\FFT_functions";

fs = 1000;  % 샘플링 주파수
N = 4096;   % 신호 길이 (FFT 포인트 수)
t = (0:N-1)/fs;

% 이상적인 사인파 신호 생성
f0 = 50; % 신호 주파수
signal = sin(2 * pi * f0 * t).*10;

% Blackman 윈도우 생성 및 적용
window = blackman(N)';
windowed_signal = signal .* window;

% FFT 계산 및 정규화
FFT_signal = fft(signal);
FFT_windowed_signal = fft(windowed_signal);
NFFT = length(FFT_signal);

% 주파수 BIN 넓이
f_bin_width = fs / NFFT;

% Power Spectrum 계산
PowerSpectrum = (abs(FFT_windowed_signal ./ N .* 2).^2);

% 주파수 BIN 설정 (f0 근처 주파수 성분 찾기)
f_bin = round(f0 / f_bin_width);

% FFT 값의 합산 (절대값)
fft_sum = sum(abs(FFT_windowed_signal(f_bin-5:f_bin+5))) .* 2 / NFFT;

% 새로운 윈도우 에너지 보정 계수 계산
window_energy_correction = fft_weight(window, NFFT);

% 주파수 BIN 합산 (넓게 설정, 여기서는 ±5 BIN)
power_sum = sum(PowerSpectrum(f_bin-5:f_bin+5));
power_sqrt = sqrt(power_sum);

% 보정된 Power Spectrum의 제곱근 계산
corrected_power_sqrt = power_sqrt.*(window_energy_correction);

% 결과 출력
printf('Frequency BIN width: %.3f Hz\n', f_bin_width);
printf('Square root of the sum of Power Spectrum: %.3f\n', power_sqrt);
printf('Corrected square root of the sum of Power Spectrum: %.3f\n', corrected_power_sqrt);
printf('Sum of the absolute FFT values: %.3f\n', fft_sum);