Matlab中的频率计算和快速傅里叶变换

Question

我有一个与快速傅里叶变换相关的问题。我想计算相位并使 FFT 绘制功率谱密度。但是，当我计算频率 f 时，出现了一些错误。这是我的程序代码：

n = 1:32768;

T = 0.2*10^-9; % Sampling period

Fs = 1/T; % Sampling frequency

Fn = Fs/2; % Nyquist frequency

omega = 2*pi*200*10^6; % Carrier frequency

L = 32768; % % Length of signal

t = (0:L-1)*T; % Time vector

x_signal(n) = cos(omega*T*n + 0.1*randn(size(n))); % Additive phase noise (random)

y_signal(n) = sin(omega*T*n + 0.1*randn(size(n))); % Additive phase noise (random)

theta(n) = atan(y_signal(n)/x_signal(n));

f = (theta(n)-theta(n-1))/(2*pi)

Y = fft(f,t);

PSD = Y.*conj(Y); % Power Spectral Density

%Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector

Answer 1

正如发布的那样，您会收到错误

error: subscript indices must be either positive integers less than 2^31 or logicals

指的是操作 theta(n-1) 当 n=1 导致索引为 0（这是越界，因为 Matlab 使用基于 1 的索引）。为避免这种情况，可以使用 n:

中的索引子集

f = (theta(n(2:end))-theta(n(1:end-1)))/(2*pi);

就是说，如果您这样做是为了尝试获得频率的瞬时测量值，那么您将有更多的问题需要处理。最微不足道的是你还应该除以 T。由于使用了 / 运算符，theta 是一个标量（参见 Matlab's mrdivide) rather than the ./ operator 执行逐元素除法。因此，更好的表达式是:

theta(n) = atan(y_signal(n)./x_signal(n));

现在，您可能会注意到的下一个问题是您实际上丢失了一些相位信息，因为 atan 的结果是 [-pi/2,pi/2] 而不是完整的 [-pi,pi] 范围。为避免这种情况，您应该改用 atan2:

theta(n) = atan2(y_signal(n), x_signal(n));

即便如此，您可能会注意到，只要相位在 -pi 附近和 pi 附近跳跃，估计的频率就会定期出现尖峰。这可以通过计算相位差模 2*pi:

来避免

f = mod(theta(n(2:end))-theta(n(1:end-1)),2*pi)/(2*pi*T);

最后要注意的一点是：调用 fft 时，您不应传入时间变量（默认假定输入以固定时间间隔进行采样）。但是，您可以指定所需的 FFT 长度。因此，您将计算 Y 如下：

Y = fft(f, L);

然后您可以使用以下方法绘制结果 PSD：

Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector
plot(Fv, abs(PSD(1:L/2+1)));

Matlab中的频率计算和快速傅里叶变换

Frequency computation and fast fourier transform in Matlab

for-loop

fft