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

当 n=1 导致索引为 0 时引用操作 theta(n-1)（这是超出范围的，因为 Matlab 使用基于 1 的索引）。为了避免这种情况，可以使用n 中的索引子集：

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

也就是说，如果您这样做是为了尝试获得频率的瞬时测量值，那么您将需要处理更多问题。最简单的一个是你也应该除以T。不那么明显的事实是，theta 是一个标量，因为使用了/ 运算符（请参阅Matlab's mrdivide），而不是执行逐元素除法的./ 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)));