2

我创建了一个与matlab 帮助页面上完全相同的对数啁啾。

t = 0:0.001:10;      % 10 seconds @ 1kHz sample rate
fo = 10; f1 = 400;   % Start at 10Hz, go up to 400Hz
X = chirp(t,fo,10,f1,'logarithmic');
figure(2);
spectrogram(X,256,200,256,1000,'yaxis');

频谱图

然后我使用以下代码将它带到频域，该代码适用于我的其他应用程序。

fft_prep = fftshift(fft(X));
fft_mag = abs(fft_prep);
pos_fft = fft_mag(1:ceil(length(fft_mag)/2));
db_fft = 20*log10(pos_fft);
figure(1);
plot(db_fft);

我很惊讶地看到下图在 1kHz-5kHz 时似乎令人兴奋：

频谱分析仪

我对 matlab 中的啁啾函数不太熟悉，想知道是否有人看到我遗漏的明显内容。欢迎任何其他指针。

4

4 回答 4

7

功能没什么问题chirp...

您只需要根据频率值绘制 db_fft，而不是矢量索引 =)。

plot(linspace(fo,f1,length(db_fft)), db_fft);

在此处输入图像描述

我还测试了使用我的其他 FFT 方法计算信号的 FFT，它们也指示了 0 到 400 Hz 之间的范围。

更新：

IMO，我发现不以 dB 或功率（周期图）绘制在视觉上更容易。这是一个很好的例子，也是我计算时域信号 FFT 的 goto 方法：mathworks.se/help/matlab/ref/fft.html

回复：

经过一番思考，我同意我上面的回答是不正确的，但不是因为你说的原因。频域中的 x 轴不应达到啁啾的实际长度（或一半，或 dubbel 或类似的东西）。频域中的 x 轴应该达到信号采样率 (Fs/2) 的一半，并且您有义务确保您以两倍于您希望/希望的最大频率的采样频率对信号进行采样解决。

换句话说，假设您的 FFT 与您的时域信号的长度相同/两倍/一半是不正确的，因为我们可以选择任意数量的频率区间来表示 FFT，最佳实践是长度 = N^2 ( 2) 的幂，用于快速计算。想一想，为什么在计算 FFT 时还需要知道时间值？你没有！您只需要采样频率（应设置为 Fs = 1000 btw，而不是 Fs = 0.001）。

我上面的回答是不正确的，应该是：

plot(linspace(0, Fs/2, length(db_fft)), db_fft)

而不是 Fs/2，你写的是 length(t)/(2*Tfinal)。它（几乎）与 Fs/2 的值相同，但它不是正确的方法 =)。

这是我的 goto FFT 方法（值不是以 dB 为单位）。

function [X,f] = myfft(x,Fs,norm)
    % usage: [X, f] = myfft(x,Fs,norm);
    %        figure(); plot(f,X)
    % norm: 'true' normalizes max(amplitude(fft))=1, default=false.
    if nargin==2
        norm=false;
    end
    L = length(x); NFFT = 2^nextpow2(L);
    f = Fs/2*linspace(0,1,NFFT/2+1);
    %f =0:(Fs/NFFT):Fs/2;
    X = fft(x,NFFT)/L; X = 2*abs(X(1:NFFT/2+1));
    if norm==true; X = X/max(abs(X)); end
end

这是 [Xfft, f] = myfft(X,Fs); 的结果图 情节（f，Xfft）；请注意，根据 NyQuist 定理，返回频率 bin 向量具有 max(f) = Fs/2（无法解析任何高于 Fs/2 的频率）。

在此处输入图像描述

于 2013-03-14T07:22:54.217 回答
0

我有几个错误，但并非所有错误都已修复。在尝试了更多代码之后，这是我想出的解决方案。

我添加了更多变量以使设置更加模块化。

Tfinal = 10;
Fs = 1/1000;
t = 0:Fs:Tfinal;      % 10 seconds @ 1kHz sample rate
fo = 10; f1 = 400;   % Start at 10Hz, go up to 400Hz
X = chirp(t,fo,Tfinal,f1,'linear');

当我绘制幅度与 linspace 的关系时，linspace 需要匹配实际啁啾信号的长度，而不仅仅是从低频到高频。因为向量t的长度为 1000，并且在 FFT 之后我们使用正半部分，所以啁啾信号的长度为 500，而不是 linspacef1建议的 400。

fft_prep = fftshift(fft(X));
fft_mag = abs(fft_prep);
pos_fft = fft_mag(ceil(length(fft_mag)/2)+1:length(fft_mag));
db_fft=20*log10(pos_fft);
figure(1);
plot(linspace(0,length(t)/(2*Tfinal),length(db_fft)), db_fft);

我还进行了前面提到的修复，以获得 FFT 的正面而不是负面。这情节：

在此处输入图像描述

这准确地描绘了 10Hz-400Hz 的啁啾激励。通过极端情况并使其线性化，您可以更清楚地看到它。尝试采样频率为 100，范围为 10-25，不进行 linspace 校正：

在此处输入图像描述

修正后：

在此处输入图像描述

于 2013-03-15T19:02:50.683 回答
0

顺便说一句，从 FFT 中提取原始啁啾幅度可能很有用。例如，如果您正在对某个设备的频率进行音频扫描，并且您想知道每个频率处响应的幅度和相位（就像 Pspice 给您的那样）。简而言之，幅度是 a^2 = abs(FFT)^2 *4 * 线性调频带宽 /(Fs * N) 其中 Fs 是采样频率，N 是 FFT 中的点数。例如，从 200 到 400Hz 的啁啾的带宽是 200Hz。

如果你想知道导数，从 Parseval 的定理开始：时间信号的均值 sqr = PSD 下的面积。因此，a^2/2 = sum(abs(FFT)^2) / N^2 其中a 是扫频信号（例如啁啾）的峰值幅度。如果啁啾在频率分布上是线性的而不是对数的，那么 FFT 是平坦的，如上图所示。因此，我们可以用 Nb * abs(FFT)^2 / N^2 替换总和，其中 Nb 是啁啾占据的频率区间数，abs(FFT) 是 FFT 的幅度，对于所有啁啾占据的频率区间。利用一个频点的带宽为 Fs/N 这一事实，我们得到啁啾的带宽为 Nb * Fs /N。上面的结果现在很容易从这里推导出来。

于 2014-05-16T02:58:49.340 回答
0

MATLAB 的文档实际上提供了fft简单的指令，这些指令通常适用于任何选择chirp（例如，二次或线性）：

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

要分析的信号示例：

S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));

计算 FFT：

Y   = fft(X);       % Calculate FFT
P2  = abs(Y/L);
P1  = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

频率向量计算为

f = Fs*(0:(L/2))/L;

如果您想生成波德幅度图，您首先应将频率转换为 [rad/s]，并将 fft 的结果转换为 [dB]：

fRad = f*2*pi;
Pdb = 20*log10(P1);

然后制作一个波特图（考虑到结果的潜在噪声性质，我建议使用散点图）

figure
scatter(fRad,Pdb)
set(gca,'xscale','log') 
grid on; grid minor
xlabel('Frequency [rad/s]')
ylabel('Magnitude [dB]')
于 2017-09-20T10:27:17.633 回答