FFT of a frequency sweep using logarithmic spacing.

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

[Pat Gipper]

Over the total duration of the chirp signal, the samples are only within a given frequency band for a small duration of the total sweep period. With the linear sweep you see an amplitude of 0.008 instead of 0.5 because the sweep occupied a given frequency band for only 1.6% (=100*0.008/0.5), whereas it has a zero amplitude value for the remaininder of the sweep period. I'm sure you can come up with the correct normalization factor with this understanding.

You have exactly the same problem with the logarithmic sweep, but the normalization factor is a bit more complicated. Good luck!

[Nuno Rocha]

Hello Pat, thanks for providing the correct answer for the linear case!

I'll keep the post open for the lograithmic spacing question, but will flag your previous reply.

Thanks,

Nuno

[Pat Gipper]

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]

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.

Answer 1

[Pat Gipper]

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)|');

[Nuno Rocha]

Hey Pat,

This is exactly what I was looking for. It works perfectly.

Thank you very much!

Nuno

[Pat Gipper]

Your welcome Nuno. I really wanted to give a complete answer to your original question, and not leave things half done.

Answer 2

[Pat Gipper]

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)|');

[Rena Berman]

@Nuno Rocha, I moved the comment to an answer. Please feel free to accept it.

[Nuno Rocha]

@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]

@Nuno Rocha Sorry no, but you can vote for the answer.