추천 게시물
FFT Peak 주파수 추정 방식 (이차 보간법, 가중평균법, parabolic inpterpolation, weighted average) 비교
- 공유 링크 만들기
- 이메일
- 기타 앱
목차
FFT를 그냥 진행하면 Leakage(누설)가 발생한다. 누설이라는 용어보단 Leakage가 편해서 Leakage로 표현하겠다. 이는 우리가 흔히 FBIN이라고 하는 주파수 Resolution이 입력 신호와 맞지 않아 발생하는 흔한 현상이다. 이해를 돕기 위해서 부연하면 FBIN은 FFT의 X축의 간격이다.
이런 현상을 완화시키기 위해서 윈도우를 씌우기도 하지만 Leakage를 완전히 막지는 못한다. 따라서 정확한 Peak 주파수를 얻기 위해선 Peak 주파수 주변 신호를 활용해 정확도를 높이는 방식을 적용해 볼 수 있다.
여기서 소개할 방식은 2가지로 신호에서 흔히 사용되는 이차 보간법과 가중평균법이다. 결론적으로 주파수를 가변하며 FFT를 하고 FFT의 결과에서 2가지 방식을 통해 주파수를 다시 계산했다. 그리고 원래 입력 주파수와 계산된 주파수의 차이를 Error로 하여 결과에 Plot하는 코드를 작성했다.
코드는 OCTAVE로 작성되어서 MATLAB에서도 무난히 동작할 것이다.
코드 예제
이차보간법 코드 (Parabolic Interpolation)
% Load the signal package
pkg load signal;
% Octave Script for Precise Peak Frequency Estimation Considering Windowing Effect
% Parameters
Fs = 128; % Sampling frequency
NFFT = 2048; % FFT points
freq_range = 0.5:0.0001:4; % Frequency range to test
errors = zeros(size(freq_range)); % Initialize error array
% Window function (Hann window)
window = hann(NFFT);
for i = 1:length(freq_range)
% Generate the signal with frequency freq_range(i)
f = freq_range(i);
t = (0:NFFT-1) / Fs;
signal = sin(2 * pi * f * t);
% Apply window
windowed_signal = signal .* window';
% FFT
X = abs(fft(windowed_signal, NFFT));
X = X(1:NFFT/2+1); % Only positive frequencies
% Find peak index
[~, k_peak] = max(X);
% Quadratic interpolation to find a more precise peak location
if k_peak > 1 && k_peak < NFFT/2
alpha = X(k_peak - 1);
beta = X(k_peak);
gamma = X(k_peak + 1);
% Parabolic interpolation formula to refine peak index
delta = 0.5 * (alpha - gamma) / (alpha - 2 * beta + gamma);
k_interp = k_peak + delta;
else
k_interp = k_peak; % No interpolation possible at boundaries
end
% Convert to frequency
f_interp = (k_interp - 1) * Fs / NFFT;
% Calculate error
errors(i) = abs(f - f_interp);
end
% Plotting the results
figure;
plot(freq_range, errors, 'b', 'LineWidth', 2);
xlabel('Input Frequency (Hz)');
ylabel('Frequency Error (Hz)');
title('Frequency Estimation Error vs Input Frequency');
grid on;
가중 평균법 코드 (Weighted Average)
% Load the signal package
pkg load signal;
% Octave Script for Precise Peak Frequency Estimation Using Amplitude Weights
% Parameters
Fs = 128; % Sampling frequency
NFFT = 2048; % FFT points
freq_range = 0.5:0.0001:4; % Frequency range to test
errors = zeros(size(freq_range)); % Initialize error array
% Window function (Hann window)
window = hann(NFFT);
for i = 1:length(freq_range)
% Generate the signal with frequency freq_range(i)
f = freq_range(i);
t = (0:NFFT-1) / Fs;
signal = sin(2 * pi * f * t);
% Apply window
windowed_signal = signal .* window';
% FFT
X = abs(fft(windowed_signal, NFFT));
X = X(1:NFFT/2+1); % Only positive frequencies
% Find peak index
[~, k_peak] = max(X);
% Weighted average using amplitude squared
if k_peak > 1 && k_peak < NFFT/2
Fm1 = (k_peak - 2) * Fs / NFFT;
F0 = (k_peak - 1) * Fs / NFFT;
Fp1 = k_peak * Fs / NFFT;
Am1 = X(k_peak - 1);
A0 = X(k_peak);
Ap1 = X(k_peak + 1);
% Weighted frequency estimation
f_estimated = (Fm1 * Am1^2 + F0 * A0^2 + Fp1 * Ap1^2) / (Am1^2 + A0^2 + Ap1^2);
else
% Edge case, use peak bin frequency
f_estimated = (k_peak - 1) * Fs / NFFT;
end
% Calculate error
errors(i) = abs(f - f_estimated);
end
% Plotting the results
figure;
plot(freq_range, errors, 'b', 'LineWidth', 2);
xlabel('Input Frequency (Hz)');
ylabel('Frequency Error (Hz)');
title('Frequency Estimation Error vs Input Frequency Using Amplitude Weights');
grid on;
코드 실행 결과
가중 평균 방식 (Weight Average)
비교 결과
- 공유 링크 만들기
- 이메일
- 기타 앱
댓글
댓글 쓰기