FFT of a frequency sweep using logarithmic spacing.

47 次查看（过去 30 天）
Nuno Rocha
Nuno Rocha 2024-5-13
评论： Rena Berman 2024-6-6
Ran in:
Hello,
I've been trying to obtain the FFT of a frequnecy sweep performed using a logarithmic progression. The signal was generated using a waveform generator, but is similar to that obtained using the chirp function as in the example below.
% Define parameters
Fs = 10000; % Sampling frequency (Hz)
T = 1/Fs; % Sampling period
L = Fs*2; % Length of signal
t = (0:L-1)*T; % Time vector
% Generate frequency sweep signal
f1 = 1; % Start frequency (Hz)
f2 = 2000; % End frequency (Hz)
y = 0.5*chirp(t, f1, T*L, f2, 'logarithmic');
% Calculate FFT
Y = fft(y);
% Calculate the frequency axis
f = Fs*(0:(L/2))/L;
% Normalize the FFT result by the amplitude of the original signal
P2 = abs(Y/L);
P1_arb = P2(1:L/2+1);
P1_arb(2:end-1) = 2*P1_arb(2:end-1);
% Plot signal and FFT
figure;
subplot(1, 2, 1)
plot(t, y)
title('Generated Signal')
xlabel('Time (s)');
ylabel('Amplitude (V)');
subplot(1, 2, 2)
plot(f, P1_arb);
title('Log Chirp');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');
When plotting the signal it's clear that the amplitude is 0.5. However, the FFT shows a variable amplitude which is nowhere close to the value of 0.5 expected.
If instead the chirp function is set to 'linear' the result is a constant amplitude across the frequency range but the amplitude is 0.008, which I don't understand how it is related to the inital value of 0.5.
y = 0.5*chirp(t, f1, T*L, f2, 'linear');
% Calculate FFT
Y = fft(y);
% Calculate the frequency axis
f = Fs*(0:(L/2))/L;
% Normalize the FFT result by the amplitude of the original signal
P2 = abs(Y/L);
P1_arb = P2(1:L/2+1);
P1_arb(2:end-1) = 2*P1_arb(2:end-1);
% Plot signal and FFT
figure;
subplot(1, 2, 1)
plot(t, y)
title('Generated Signal')
xlabel('Time (s)');
ylabel('Amplitude (V)');
subplot(1, 2, 2)
plot(f, P1_arb);
title('Linear Chirp');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');
Could you help me normalize the resulting FFT in order to obtain the expected amplitude.
Thank you,
Nuno
  4 个评论
显示 2更早的评论
Pat Gipper
Pat Gipper 2024-5-14
Thanks Nuno. I found the solution to the linear case using emperical methods and do not have an analytical means to get to the answer that I did. Using the same method for the logarithmic case would involve a lot of trial and error. So I'll leave that one for some smart person out there!
Nuno Rocha
Nuno Rocha 2024-5-15
Thanks for your help Pat. That was a very clever find nonetheless!
I'll post the normalization method for the logarithmic case here if I eventually find it.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Pat Gipper
Pat Gipper 2024-5-17
Here is some code to perform the logarithmic chirp normalization.
% Define parameters
clear norm
N = 1; % Chirp duration (sec)
Fs = 20000; % Sampling frequency (Hz)
T = 1/Fs; % Sampling period
L = Fs*N; % Length of signal
t = (0:L-1)*T; % Time vector
% Generate frequency sweep signal
f1 = 1; % Start frequency (Hz)
f2 = 2000; % End frequency (Hz)
y = 0.5*chirp(t, f1, T*L, f2, 'logarithmic');
% Calculate FFT
Y = fft(y);
% Calculate the frequency axis
f = Fs*(0:(L/2))/L;
% This section of code calculates a correction factor for a given fft bin
% whereby it adjusts inversely for the time the chirp is within that bin
% frequency range versus the total chirp duration
B = (f2/f1)^(1/N); % chirp constant for logarithmic sweep
for i = 1:size(f,2)-1
flow = f(i);
fhigh = f(i+1);
thigh = log((fhigh/f1))/log(B);% Solve for the end time of this fft bin
tlow = log((flow/f1))/log(B);% Solve for the start time of this fft bin
norm(i) = N / (thigh - tlow);% Normalization factor of this fft bin
end
% Normalize the FFT result by the amplitude of the original signal
P2 = abs(Y/L);
P1_arb = P2(1:L/2+1);
P1_arb(2:end-1) = 2*P1_arb(2:end-1);% Adjust non-DC terms for negative freq
P1_arb(1:end-1) = P1_arb(1:end-1) .* sqrt(norm);% Logrithmic normalization
% Plot signal and FFT
figure;
subplot(1, 2, 1)
plot(t, y)
title('Generated Signal')
xlabel('Time (s)');
ylabel('Amplitude (V)');
subplot(1, 2, 2)
plot(f, P1_arb);
title('Log Chirp');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');
  2 个评论
Nuno Rocha
Nuno Rocha 2024-5-20
Hey Pat,
This is exactly what I was looking for. It works perfectly.
Thank you very much!
Nuno
Pat Gipper
Pat Gipper 2024-5-20
Your welcome Nuno. I really wanted to give a complete answer to your original question, and not leave things half done.

请先登录，再进行评论。

更多回答（1 个）

Pat Gipper
Pat Gipper 2024-5-14
移动：Rena Berman 2024-6-5
The correction factor for a linear chirp is sqrt(beta) / df, where beta is the chirp rate in Hz/sec and df is the FFT bin size in Hz. Here is the modified code that includes this added correction factor.
% Define parameters
N = 2; % Chirp duration (sec)
Fs = 10000; % Sampling frequency (Hz)
T = 1/Fs; % Sampling period
%L = Fs*2; % Length of signal
L = Fs*N; % Length of signal
t = (0:L-1)*T; % Time vector
% Generate frequency sweep signal
f0 = 1; % Start frequency (Hz)
f1 = 2000; % End frequency (Hz)
y = 0.5*chirp(t, f0, T*L, f1, 'linear');
% What is the sweep rate Beta
beta = (f1-f0)/max(t); % f(t) f1+beta*t
% Calculate FFT
Y = fft(y);
% Calculate the frequency axis
f = Fs*(0:(L/2))/L; % Only plot out to the Nyquist rate
df = f(2); % Frequency bin width (Hz)
% Normalize the FFT result by the amplitude of the original signal
P2 = abs(Y/L); % Normalize by the number of samples
P1_arb = P2(1:L/2+1); % Plot out to the Nyquist rate
P1_arb(2:end-1) = 2*P1_arb(2:end-1); % Double to include negative freqs
P1_arb(2:end-1) = P1_arb(2:end-1) * sqrt(beta) / df; % Correction factor
% Plot signal and FFT
figure;
subplot(1, 2, 1)
plot(t, y)
title('Generated Signal')
xlabel('Time (s)');
ylabel('Amplitude (V)');
subplot(1, 2, 2)
plot(f, P1_arb);
title('Linear Chirp');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');
  3 个评论
显示 1更早的评论
Nuno Rocha
Nuno Rocha 2024-6-6
@Rena Berman, thanks for doing this. In the meantime I've accepted another answer by the same user which answers the main question posted. This reply is to a secondary question from my original post. Is there a way to accept both?
Thank you!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by