Matlab绘制啁啾的对数幅度答案

【问题标题】：Matlab plot log magnitude of chirpMatlab绘制啁啾的对数幅度
【发布时间】：2013-03-14 04:18:31
【问题描述】：

我创建了一个与matlab help page 完全相同的对数啁啾。

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 中的啁啾函数不太熟悉，想知道是否有人看到我遗漏的明显内容。欢迎任何其他指针。

【问题讨论】：

对啁啾信号进行 FFT 似乎不是一件有用的事情？我不认为只是因为您的频谱图的范围是 10Hz-400Hz，而这些频率与波特图上显示的频率相同。我认为您的问题更多地与有限时间啁啾信号的频率特性有关，而不是与 matlab 有关。我建议你在这里问：dsp.stackexchange.com
好点，也许那个频谱图没有向我展示整个画面。这是我试图更好地理解的啁啾函数。我希望我可以扫描有限数量的频率并通过 FFT 看到它们。你认为我希望看到什么？我使这个问题更加以代码为中心，专注于可能在我的啁啾函数或其他东西中找到错误。我也在 DSP 论坛上发布了更多 DSP-based question。如果您认为它可以更友好，请随时编辑我的问题。
不，我认为它不属于这里，您的问题与代码无关。啁啾信号在不同时间显示不同的频率，但 FFT 看起来是整个信号（所有时间），为了实现这一点，我想它会产生更高的谐波。 matlab 示例绘制频谱图并且不采用 FFT 是有原因的，这是因为啁啾的 FFT 并不是真正有用的事情。啁啾信号本身在一定频率范围内测试系统，而不必离开时域，重点（在我看来）是消除进行 FFT 的需要
顺便说一句，如果您要在有限的频率范围内寻找恒定幅度（即频域中的 rect 函数），则转换后的时域函数是 sinc 函数，即 @ 987654330@.
我绝对同意你的观点，对啁啾进行 FFT 并不是一件很有用的事情。我正在尝试详细了解啁啾功能，并想验证我正在用它创建什么。我试图将这个理论排除在这个线程之外，但我们的 cmets 肯定让我们误入歧途，也许我们应该将我们的 cmets 移到 DSP 线程或将它们带到电子邮件通信中。正如下面的答案所证实的那样，我的代码中有一个错误或误解，我试图关注这些问题，我觉得这个线程对我来说非常有价值，希望将来对其他人也很有价值。

标签： matlab signal-processing fft frequency-analysis

【解决方案1】：

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)

你写的是 length(t)/(2*Tfinal) 而不是 Fs/2。它（几乎）与 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 的频率都无法解析）。

【讨论】：

我可能弄错了，但 fft 信号不是仍以左右镜像方式绘制吗？使用fftshift 后，频谱以零为中心，当您仅绘制左半部分时，您的绘图右侧最终为零，对吧？
确实，我也注意到了。
我假设您在谈论这一行：pos_fft = fft_mag(1:ceil(length(fft_mag)/2)); 切换到pos_fft = fft_mag(ceil(length(fft_mag)/2)+1:length(fft_mag));？
是的，或者您可以简单地使用 fft(fftshift(X))，而不是 fftshift(fft(X))
感谢您帮助我找到这些编码问题！我已经更新了这个DSP thread 来询问更多关于理论的信息并学习如何创建一个恒定幅度的啁啾声。

【解决方案2】：

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

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

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 而不是 400，直到 f1 的 linspace 会建议。

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 校正：

修改后：

【讨论】：

【解决方案3】：

顺便说一下，从 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。上面的结果现在很容易从这里推导出来。

【讨论】：

【解决方案4】：

MATLAB 的 documentation 关于 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]')

【讨论】：