FFT 예제: 윈도우 적용 시 Normalized Wegith Value, SNR 계산, Sinusoidal + Noise Signal Generation, Prime Cycle Frequency 등 포함

목차

측정한 신호를 FFT로 변환하기 위해서는 다양한 방법들을 고려 해야 한다. 단순히 인터넷에 있는 ideal sinewave를 바탕으로 깔끔하게 나온 FFT 결과는 현실적으로 적용하기 어렵기 때문이다.

아래 항목들이 측정 및 FFT 코드 작성 시 고민해야할 사항이다.

1. Normalize

2. Sine wave cycle

3. Signal and Noise power, rms calculation

4. Window selection and weight

각각의 항목에 관해서는 나중에 별도의 글로 다시 정리할 예정이다.

아래는 위의 사항을 고려하여 작성한 FFT 코드이다.

OCTAVE를 기준으로 작성되었다.

FFT 코드

clc;

pkg load signal;

pkg load io;

% Functions

function weight_result = fft_weight2(window, NFFT)

  %weight_result=1/sum(abs(fft(window,NFFT)./NFFT.*2));

  window_fft=abs(fft(window,NFFT)./NFFT);

  weight_result=1/sqrt(sum(window_fft.^2))

  %window2=window.^2;

  %weight_result=sqrt(sum(abs(fft(window2,NFFT)./NFFT)));

  return;

end

function adjusted_frequency = closest_prime_frequency(sampling_time, target_frequency)

  % Calculate the number of cycles for the target frequency

  num_cycles = sampling_time * target_frequency;

  % Define the search range around num_cycles

  search_range = 1000;

  % Calculate the lower and upper bounds for prime search

  lower_bound = max(2, floor(num_cycles - search_range));

  upper_bound = ceil(num_cycles + search_range);

  % Generate the list of prime numbers within the search range

  primes_list = primes(upper_bound);

  primes_list = primes_list(primes_list >= lower_bound);

  % Find the prime number closest to the calculated number of cycles

  [~, idx] = min(abs(primes_list - num_cycles));

  closest_prime = primes_list(idx);

  % Calculate the adjusted frequency

  adjusted_frequency = closest_prime / sampling_time;

  % Display the results

  %{

  fprintf('Original Frequency: %.8f Hz\n', target_frequency);

  fprintf('Sampling Time: %.8f seconds\n', sampling_time);

  fprintf('Closest Prime Number of Cycles: %d\n', closest_prime);

  fprintf('Adjusted Frequency: %.8f Hz\n', adjusted_frequency);

%}

  % Return the adjusted frequency

  return;

end

# Functions Ends

% Parameters

Fs = 256;                 % 샘플링 주파수

T = 1/Fs;                 % 샘플링 주기

NFFT = 2^14;                 % 신호 길이 (NFFT)

t = (0:NFFT-1)*T;            % 시간 벡터

hanning_window =  hann(NFFT)'; % weight 1.633

hamming_window = hamming(NFFT)'; % weight 1.58635

flattop_window = flattopwin(NFFT)'; % weight 2.3898

blackman_window = blackman(NFFT)'; % wieght 1.812

%time domain window weight calculation

hanning_weight=1/((sum(hanning_window.^2)./NFFT)^0.5)

hamming_weight=1/((sum(hamming_window.^2)./NFFT)^0.5)

flattop_weight=1/((sum(flattop_window.^2)./NFFT)^0.5)

blackman_weight=1/((sum(blackman_window.^2)./NFFT)^0.5)

%frequency domain window wieght calculation

hanning_weight2=fft_weight2(hanning_window,NFFT)

hamming_weight2=fft_weight2(hamming_window,NFFT)

flattop_weight2=fft_weight2(flattop_window,NFFT)

blackman_weight2=fft_weight2(blackman_window,NFFT)

sig_rms_arr=zeros(2,100);

noise_rms_arr=zeros(2,100);

for i=1:1 % Freq.에 따라 window 적용 후 power 일정함 확인

% Signal A V @ f0 Hz sinewave

f0 = i;                   % 신호 주파수

A = 10;                % 10V peak 진폭

fsig=closest_prime_frequency(T*NFFT, f0);

signal = (A) * sin(2*pi*fsig*t);

% Noise B: noise_rms Vrms 랜덤 노이즈

noise_rms = 1;        % 1Vrms

noise = noise_rms * randn(size(t));

noise_std = std(noise);

% Combined signal

combined_signal = signal + noise;

combined_signal_win = combined_signal .* hanning_window.*1.633;

%combined_signal_win = combined_signal .* hamming_window.*1.5859;

%combined_signal_win = combined_signal .* flattop_window.*2.387;

%combined_signal_win = combined_signal .* blackman_window.*1.81;

% FFT

%Y = fft(combined_signal, NFFT);

Y = fft(combined_signal_win, NFFT);

P2 = abs(Y/NFFT);            % 양면 스펙트럼

P1 = P2(1:NFFT/2);         % 단일면 스펙트럼

P1(2:end-1) = P1(2:end-1).*sqrt(2);

% 주파수 벡터

f = Fs*(0:(NFFT/2)-1)/NFFT;

% 신호 파워 계산

signal_fbin = 1;         % FBIN = ±20

signal_power = sum(P1((f >= fsig-signal_fbin) & (f <= fsig+signal_fbin)).^2);

% 노이즈 파워 계산을 위해 신호 영역을 0으로 설정

P1_no_signal = P1;

P1_no_signal((f >= fsig-(signal_fbin*Fs/NFFT )) & (f <= fsig+(signal_fbin*Fs/NFFT))) = 0;

% 노이즈 파워 계산

noise_power = sum(P1_no_signal.^2);

% SNR 계산

SNR = 10 * log10(signal_power / noise_power);

sig_rms_arr(1,i)=fsig;

sig_rms_arr(2,i)=sqrt(signal_power);

noise_rms_arr(1,i)=fsig;

noise_rms_arr(2,i)=sqrt(noise_power);

 if i==1

 % 결과 출력

 fprintf('Signal Power: %.4e W\n', signal_power);

 fprintf('Signal Rms: %.4e W\n', sqrt(signal_power));

 fprintf('Noise Power: %.4e W\n', noise_power);

 fprintf('Noise Rms: %.4e W / STDEV(time domain): %.4e W\n', sqrt(noise_power), noise_std);

 fprintf('SNR: %.2f dB\n', SNR);

 % 신호의 시간축 그래프

 figure(1);

 plot(t, combined_signal_win);

 title('Time-Domain Signal');

 xlabel('Time (s)');

 ylabel('Amplitude (V)');

 % FFT 그래프 (로그 스케일)

 figure(2);

 semilogx(f(2:(NFFT/2)), 20*log10(P1(2:(NFFT/2))),'o-');

 title('Frequency-Domain Signal (Log Scale)');

 xlabel('Frequency (Hz)');

 ylabel('Amplitude (dB)');

 xlim([0, Fs/2]);

 % 그래프 표시

 figure(gcf);

 end

end

코드 간략 설명

window의 weight 값을 time domain과 frequency domain에서 확인. 

입력 주파수를 소수 사이클을 갖는 가장 가까운 주파수로 변경.

입력 신호와 노이즈를 생성 (값은 설정 가능).

noise(time domain)의 Standard deviation 계산.

신호 + Noise 결합

FFT 후 Nyquist 주파수 영역 만 취한 뒤 Normalize

fbin range를 설정하고 해당 영역 신호의 power 계산

fbin 안쪽의 신호 영역을 0으로 한 뒤 나머지 영역 합산하여 노이즈 계산

파형 plot

예제에서 신호의 Amplitude는 10V이므로, 10 * sine(t)를 제곱하면 100 * (1/2 - cos(2t)/2)가 된다.

0~2pi까지 적분한 뒤 2pi 로 나누면 평균 값이 된다.

100 * t/2 - 100/4 * sine(2t) + C

t에 0과 2pi를 넣고 빼준다. sine과 C는 사라진다.

100 * pi

적분 구간으로 나눠 평균값을 구한다.

100 * pi / (2 pi) = 50

power는 50이 되며, RMS는 7.07이 된다. 즉 rms는 amplitude / sqrt(2) 와 같다.