Borrar filtros
Borrar filtros

FFT of a frequency sweep using logarithmic spacing.

116 visualizaciones (últimos 30 días)
Nuno Rocha
Nuno Rocha el 13 de Mayo de 2024
Comentada: Rena Berman el 6 de Jun. de 2024 a las 12:12
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 comentarios
Pat Gipper
Pat Gipper el 14 de Mayo de 2024 a las 18:41
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 el 15 de Mayo de 2024 a las 9:16
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.

Iniciar sesión para comentar.

Respuesta aceptada

Pat Gipper
Pat Gipper el 17 de Mayo de 2024 a las 22:05
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 comentarios
Nuno Rocha
Nuno Rocha el 20 de Mayo de 2024 a las 9:16
Hey Pat,
This is exactly what I was looking for. It works perfectly.
Thank you very much!
Nuno
Pat Gipper
Pat Gipper el 20 de Mayo de 2024 a las 12:43
Your welcome Nuno. I really wanted to give a complete answer to your original question, and not leave things half done.

Iniciar sesión para comentar.

Más respuestas (1)

Pat Gipper
Pat Gipper el 14 de Mayo de 2024
Movida: Rena Berman el 5 de Jun. de 2024 a las 17:38
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 comentarios
Nuno Rocha
Nuno Rocha el 6 de Jun. de 2024 a las 8:29
@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!
Rena Berman
Rena Berman el 6 de Jun. de 2024 a las 12:12
@Nuno Rocha Sorry no, but you can vote for the answer.

Iniciar sesión para comentar.

Productos


Versión

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by