每次遇到Fourier变换都忘记的基本代码。
如何获得“右半边”的Fourier变换谱?
所谓“右半边”的Fourier变换谱其实就是频率大于零那半边。对FFT来说其实是截取“左半边”。
function [X,freq]=lhsFFT(x,Fs)
% Cuts the left-hand side of FFT
% x is the signal sample to be FTT-ed, Fs is the sampling frequency of the signal
% X is the FFT result, freq is the frequency axis
N = length(x);
freq_bin = 0:N-1; % the frequency bin
freq_res = Fs / N; % the frequency resolution
freq = freq_bin * freq_res; %create the frequency range
X = 2 * fft(x) / N; % normalize the data
cutoff = ceil(N/2); %length of the left-hand side
freq = freq(1:cutOff);
X = X(1:cutOff);
获取某频率的幅值
问题其实就是决定某频率值落在第几个frequency bin,然后四舍五入取最接近的那个频率。
N = length(x);
freq_Res = Fs / N; % the frequency resolution
[X, freq] = lhsFFT(x, Fs) % get the left-hand side FFT
f = 1.234 % The queried frequency is 1.234 Hz for example
freq_bin_prime = round(f/freq_res); % finds the frequency bin closest to f
freq_prime = freq_bin * freq_res; % the resolved frequency closest to f
X_prime = X(freq_bin_prime + 1); % the amplitude of the resolved frequency (vector index starts from 1)