在 MATLAB 中使用 FFT 的频率响应答案

【问题标题】：Frequency response using FFT in MATLAB在 MATLAB 中使用 FFT 的频率响应
【发布时间】：2010-10-22 00:56:49
【问题描述】：

这是场景：使用频谱分析仪，我有输入值和输出值。采样数为32000，采样率为2000samples/sec，输入为50 hz的正弦波，输入为电流，输出为压力，单位为psi。

我如何使用 MATLAB 从这些数据中计算频率响应，在 MATLAB 中使用 FFT 函数。

我能够生成一个正弦波，它给出幅度和相位角，这是我使用的代码：

%FFT Analysis to calculate the frequency response for the raw data
%The FFT allows you to efficiently estimate component frequencies in data from a discrete set of values sampled at a fixed rate

% Sampling frequency(Hz)
Fs = 2000;   

% Time vector of 16 second
t = 0:1/Fs:16-1;   

% Create a sine wave of 50 Hz.
x = sin(2*pi*t*50);                                                       

% Use next highest power of 2 greater than or equal to length(x) to calculate FFT.
nfft = pow2(nextpow2(length(x))) 

% Take fft, padding with zeros so that length(fftx) is equal to nfft 
fftx = fft(x,nfft); 

% Calculate the number of unique points
NumUniquePts = ceil((nfft+1)/2); 

% FFT is symmetric, throw away second half 
fftx = fftx(1:NumUniquePts); 

% Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x
mx = abs(fftx)/length(x); 

% Take the square of the magnitude of fft of x. 
mx = mx.^2; 

% Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
% The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2.

if rem(nfft, 2) % odd nfft excludes Nyquist point
  mx(2:end) = mx(2:end)*2;
else
  mx(2:end -1) = mx(2:end -1)*2;
end

% This is an evenly spaced frequency vector with NumUniquePts points. 
f = (0:NumUniquePts-1)*Fs/nfft; 

% Generate the plot, title and labels. 
subplot(211),plot(f,mx); 
title('Power Spectrum of a 50Hz Sine Wave'); 
xlabel('Frequency (Hz)'); 
ylabel('Power'); 

% returns the phase angles, in radians, for each element of complex array fftx
phase = unwrap(angle(fftx));
PHA = phase*180/pi;
subplot(212),plot(f,PHA),title('frequency response');
xlabel('Frequency (Hz)')
ylabel('Phase (Degrees)')
grid on

我在90度相位角处从相位图中获取频率响应，这是计算频率响应的正确方法吗？

如何将此响应与从分析仪获得的值进行比较？这是一个交叉检查，以查看分析器逻辑是否有意义。

【问题讨论】：

标签： matlab signal-processing fft

【解决方案1】：

乍一看还不错，但你缺少一些东西：

您应该在 FFT 之前将窗口函数应用于时域数据，例如http://en.wikipedia.org/wiki/Window_function 用于一般窗口，http://en.wikipedia.org/wiki/Hann_window 用于最常用的窗口函数（Hann aka Hanning）。
您可能希望以 dB 为单位绘制对数幅度，而不仅仅是原始幅度

【讨论】：

【解决方案2】：

您应该考虑查看用于计算频率响应的 cpsd() 函数。为您处理各种窗口函数的缩放和标准化。

频率响应将是

G = cpsd (output,input) / cpsd (input,input)

然后取angle() 获取输入和输出之间的相位差。

你的代码 sn-p 没有提到输入和输出数据集是什么。

【讨论】：