Why does a frequency input for a chirp not work?

Hello,

From references I saw that a linear sweep/chirp function creates frequencies which are not the actual frequencies which you can count if you see the output. There are formulas for calculating the right frequency for the input, but I don't understand why this is the case, as for a constant frequency this does work. See my code below for a clarification of my question

t = 0:0.01:5; % Time vector

%% Normal sine wave with 2 Hz
f=2*ones(1,length(t)); % Create a constant frequency vector
figure;
plot_y(t,f); % Plot the sine wave for this constant frequency vector

%% Linear sine wave from 1 to 2 Hz
f0 = 1; % Frequency at t=0
f1 = 4; % Frequency at t=t(end)
f=linspace(1,4,length(t)); % Create a linear frequency vector
hold on;
plot_y(t,f); % Plot the sine wave for this linear frequency vector

% This is how you should calculate a linear frequency for a sweep:
f=f0 + (f1-f0)/(2*t(end)).*t;
fprintf('f0 = %.1f, f1 = %.1f, f(1) = %.1f, f(end) = %.1f \n',f0,f1,f(1),f(end))

f0 = 1.0, f1 = 4.0, f(1) = 1.0, f(end) = 2.5 

%% Plot info
legend('normal','linear')
title('Frequency from 1 to 2 Hz and only 2 Hz')

%% Function to plot (so you can see it is indeed the same way I calculate and plot it)
function plot_y(t,f)
y = sin(2*pi.*f.*t); % Normal sine function, where the frequency can be a function of time
plot(t,y);
end

From the graph you can clearly see that the linear frequency at the end is more than 4 Hz. With the other formula where f(end) = 2.5 Hz, the function works correctly.

Could anyone explain why this is the case? The reason I need to know this is that I also want to use a polyfit to estimate a range of frequencies, which I cannot approximate with an exponential or hyperbolic fit.

Thanks in advance!

Answer 1

[Mathieu NOE]

hello

fixed a few mistakes in your code

the most obvious one is when you compute a sinewave . You have to remember that you compute the sin of an angle , whatever the frequency behaviour is (fixed or variable), the frequency is the time derivative of the angle, or vice versa the angle must be computed as the time integral of the frequency (can be time dependant as in your example)

what you did is a simple product of the time vector by the frequency (at each time index) but this is not the integral of the frequency .

the other point is why the factor 2 in this line ?

% This is how you should calculate a linear frequency for a sweep:

f=f0 + (f1-f0)/(2*t(end)).*t;

My results

Nb the computed frequency (second subplot ) matches the signal specs (from 1 to 4 Hz , linear chirp)

code fixed :

clearvars;

dt = 0.01;
t = 0:dt:5; % Time vector
%% Normal sine wave with 1 Hz
f=1*ones(1,length(t)); % Create a constant frequency vector
figure;
y = compute_y(t,f); % Plot the sine wave for this constant frequency vector
figure(1)
subplot(2,1,1),plot(t,y);

%% Linear sine wave from 1 to 4 Hz
f0 = 1; % Frequency at t=0
f1 = 4; % Frequency at t=t(end)
f=linspace(1,4,length(t)); % Create a linear frequency vector
hold on;
y = compute_y(t,f); % Plot the sine wave for this linear frequency vector

% find signal frequency by measuring "zero" (or any other value crossing
% time index

threshold = 0; % 
t0_pos1 = find_zc(t,y,threshold);
period = diff(t0_pos1); % delta time of crossing points
freq = 1./period; % signal frequency

subplot(2,1,1),plot(t,y,'r-',t0_pos1,threshold*ones(size(t0_pos1)),'.r','linewidth',2,'markersize',20);grid on
xlim([min(t) max(t)]);
legend('normal','linear sweep','linear sweep crossing points')
title('Frequency from 1 to 4 Hz and only 1 Hz')

xlabel('Time (s)');
ylabel('y');
subplot(2,1,2),plot(t0_pos1(2:end),freq,'r.-','linewidth',2,'markersize',12);grid on
xlim([min(t) max(t)]);
ylim([0 5]);
legend('signal rate (frequency)');
xlabel('Time (s)');
ylabel('Hz');

% This is how you should calculate a linear frequency for a sweep:
% f=f0 + (f1-f0)/(2*t(end)).*t; % why the factor 2 at the denominator ? 
f=f0 + (f1-f0)/(t(end)).*t; % should be this 

fprintf('f0 = %.1f, f1 = %.1f, f(1) = %.1f, f(end) = %.1f \n',f0,f1,f(1),f(end))
% f0 = 1.0, f1 = 4.0, f(1) = 1.0, f(end) = 2.5 
%% Plot info

%%% compute y (the right way)
function y = compute_y(dt,f)
    angl = cumtrapz(dt,f);
    y = sin(2*pi*angl); %
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)

    % positive slope "zero" crossing detection, using linear interpolation
    y = y - threshold;
    zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
    ix=zci(y);                      %find indices of + zero crossings of x
    ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
    Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[TG]

Thank you very much for this elaborate answer!

I used the factor 2 in the linear chirp as this is the mathematical formula for computing a sine wave without using the integral. But this is better as this will also work for a logarithmic or hyperbolic input of the frequency.

[Mathieu NOE]

Glad I could be of some help

have a great day !