Creating a sinewave with logarithmic frequency

I'm trying to generate a sinewave with logarithmic increasing frequency.

freq = f_str*((f_end/f_str).^(time/T));

The max frequency in the sinewave is 50 Hz. However, when I construct such sinewave and run a spectrogram on it (STFT), it shows that the max frequency is 150 Hz.

Any thoughts would be appreciated.

% Generate log-frequency sweep then plot it and its spectrogram

NPTS = 225000;
f_str = 7;
f_end = 50;
T = 450;

time = linspace(0,T,NPTS)'; 
freq = f_str*((f_end/f_str).^(time/T));
Y = sin(2*pi*freq.*time);       % Sinewave with logarithmic sweep
fs = (2*NPTS)/900;

subplot(311), plot(time, freq, '-b', 'LineWidth',3), grid on, xlabel('Time  [sec]'), ylabel('Frequency  [Hz]')
subplot(312), plot(time, Y, '-b', 'LineWidth',3), grid on, xlabel('Time  [sec]'), ylabel('SineWave  [-]'), axis([0 450 -2 2])
subplot(313), stft(Y, fs, 'FrequencyRange', 'onesided'), set(gca,'Ylim',[0 200]), colorbar off

[dpb]

% Generate log-frequency sweep then plot it and its spectrogram

NPTS = 225000;
f_str = 7;
f_end = 50;
T = 450;

time = linspace(0,T,NPTS)'; 
freq = f_str*((f_end/f_str).^(time/T));
Y = sin(2*pi*freq.*time);       % Sinewave with logarithmic sweep
fs = (2*NPTS)/900;

subplot(311), plot(time, freq, '-b', 'LineWidth',3), grid on, xlabel('Time  [sec]'), ylabel('Frequency  [Hz]')
subplot(312), plot(time, Y, '-b', 'LineWidth',3), grid on, xlabel('Time  [sec]'), ylabel('SineWave  [-]'), axis([0 450 -2 2])
subplot(313), stft(Y, fs, 'FrequencyRange', 'onesided'), set(gca,'Ylim',[0 200]), colorbar off
subplot(312),xlim([449.9 450])

Shows the time series generated 15 cycles/0.1 sec --> 150 Hz so stft isn't lying; you just don't have the input signal you think you have...

[Djamil Boulahbal]

That is where the puzzle is. Somehow the signal I calculated is not what it 'should' be.

The frequency (vs. time) is correct - it's the top plot.

Then I calculate the sinewave as:

Y = sin(2*pi*freq.*time); - Somehow this formula does NOT generate the correct trace

The question then is: What is so wrong with the generic expression:

sin ( 2 × pi × freq × time ) ..... ????

It seems pretty basic ... but apparently not.

[dpb]

% Generate log-frequency sweep then plot it and its spectrogram

NPTS = 225000;
f_str = 7;
f_end = 50;
T = 450;

time = linspace(0,T,NPTS)'; 
freq = f_str*((f_end/f_str).^(time/T));
Y = sin(2*pi*freq.*time);       % Sinewave with logarithmic sweep
fs = (2*NPTS)/900;
tView=[0 1]; 
y=chirp(time,f_str,T,f_end,"logarithmic",-90);
subplot(321), plot(time, Y, '-b', 'LineWidth',1), grid on, xlabel('Time  [sec]'), ylabel('Original'), axis([0 450 -2 2])
xlim(tView)
subplot(323), plot(time, y, '-b', 'LineWidth',1), grid on, xlabel('Time  [sec]'), ylabel('Chirp'), axis([0 450 -2 2])
xlim(tView)
tView=[T-0.1 T];
subplot(322), plot(time, Y, '-b', 'LineWidth',1), grid on, xlabel('Time  [sec]'), ylabel('Original'), axis([0 450 -2 2])
xlim(tView)
subplot(324), plot(time, y, '-b', 'LineWidth',1), grid on, xlabel('Time  [sec]'), ylabel('Chirp'), axis([0 450 -2 2])
xlim(tView)

Compares to using the Signal Processing TB chirp function -- I didn't have time to try to dig into just what the root cause is in the straightforward calculation; perhaps the issue of the large angle is part, but the internals of the chirp function deal with it and do produce the expected number of cycles over the last 100 msec of the time history whereas the original clearly shows that the frequency of the actual time trace itself is right at the 150 Hz, not the expected 50 Hz. OTOMH, I don't have any better explanation, only that one can illustrate the result...

[Djamil Boulahbal]

Bravo! The chirp function did the trick. I somehow had completely forgotten about it ... Thank you for your help. I am convinced now there is a numerical issue in my approach. Though no time to dwell ... as I have the answer now.

Answer 1

[Star Strider]

Using the logspace function could work.

Try this ’—

% Generate log-frequency sweep then plot it and its spectrogram

NPTS = 225000;
f_str = 7;
f_end = 50;
T = 450;

time = linspace(0,T,NPTS)'; 
% freq = f_str*((f_end/f_str).^(time/T));
freq = logspace(-3, log10(f_end/f_str),NPTS)';
Y = sin(2*pi*freq.*time);       % Sinewave with logarithmic sweep
fs = (2*NPTS)/900;

subplot(311), plot(time, freq, '-b', 'LineWidth',3), grid on, xlabel('Time  [sec]'), ylabel('Frequency  [Hz]')
subplot(312), plot(time, Y, '-b', 'LineWidth',3), grid on, xlabel('Time  [sec]'), ylabel('SineWave  [-]'), axis([0 450 -2 2])
subplot(313), stft(Y, fs, 'FrequencyRange', 'onesided'), set(gca,'Ylim',[0 200]), colorbar off

Make appropriate changes to get the result you want.

.

[Djamil Boulahbal]

Thank you. I triied it, but still it does not work. SOmehow the frequency generated in the sinewave is not as expected. In your graphs, the max frequency (top trace) is near 7 Hz. Whereas in the spectrogram (bottom colormap), the max frequency is near 70 Hz ... almost a factor of 10.

What is so wrong with the logic/formulae in here?

[Star Strider]

Looking at the various relevant arguments —

% Generate log-frequency sweep then plot it and its spectrogram

NPTS = 225000;
f_str = 7;
f_end = 50;
T = 450;

time = linspace(0,T,NPTS)'; 
% freq = f_str*((f_end/f_str).^(time/T));
freq = logspace(-3, log10(f_str),NPTS)';
Y = sin(freq.*time);       % Sinewave with logarithmic sweep
fs = (2*NPTS)/900

fs = 500

subplot(311), plot(time, freq, '-b', 'LineWidth',3), grid on, xlabel('Time  [sec]'), ylabel('Frequency  [Hz]')
subplot(312), plot(time, Y, '-b', 'LineWidth',3), grid on, xlabel('Time  [sec]'), ylabel('SineWave  [-]'), axis([0 450 -2 2])
subplot(313), stft(Y, fs, 'FrequencyRange', 'onesided'), set(gca,'Ylim',[0 20]), colorbar off

The frequency is actually correct. It goes from near D-C to 7 Hz (first plot). I looked at the stft documentatiion, and I don’t see any errors in my code or the stft call. Removing from the sin argument ‘sort of’ works.

.

[Djamil Boulahbal]

Yes, the initial frequency (7 Hz) is correct, but the frequemcy at the end of the chirp is not... at least that is what I got. That is why I believe the issue is more 'numerical' than anything else. Thanks again for the effort and answers. Matlab has a built-in chirp function, I tried it, and it works perfectly.

[Star Strider]

As always, my pleasure!

I am not cerrtain where the problem lies. The code appears to be correct, and iis by my analyses (not shown, since they did not add any useful information).

[Djamil Boulahbal]

Thanks again, you've been very thorough and helpful.

Answer 2

[Djamil Boulahbal]

Thinking a little more about it, the errror could be because the angle (2×pi×freq×time) grows very large (reaches 15×10^4 radians), and the precision is lost.

One thing I could try is to define the angle only to within +/- pi ...

Thoughts?

[dpb]

But the sin function itself deals with long time series ok it seems...

% Generate log-frequency sweep then plot it and its spectrogram

NPTS = 225000;
f_str = 7;
f_end = 50;
T = 450;

time = linspace(0,T,NPTS)'; 
%freq = f_str*((f_end/f_str).^(time/T));
%Y = sin(2*pi*freq.*time);       % Sinewave with logarithmic sweep
Y = sin(2*pi*f_str*time);        % Sinewave with starting frequency
tView=[0 1]; 
y = sin(2*pi*f_end*time);        % Sinewave with ending frequency
subplot(321), plot(time, Y, '-b', 'LineWidth',1), grid on, xlabel('Time  [sec]'), ylabel('Start F'), axis([0 450 -2 2])
tView=[0 1]; xlim(tView)
subplot(322), plot(time, y, '-b', 'LineWidth',1), grid on, xlabel('Time  [sec]'), ylabel('End F'), axis([0 450 -2 2])
tView=[T-0.1 T];xlim(tView)
subplot(323), plot(time, Y, '-b', 'LineWidth',1), grid on, xlabel('Time  [sec]'), ylabel('Start F'), axis([0 450 -2 2])
tView=[0 1]; xlim(tView)
subplot(324), plot(time, y, '-b', 'LineWidth',1), grid on, xlabel('Time  [sec]'), ylabel('End F'), axis([0 450 -2 2])
tView=[T-0.1 T];xlim(tView)

The first 1 sec shows 7 cycles, the last 100 msed is 5, so the two full-length calculations of a fixed-frequency sine wave over the same t vector produce the expected frequencies. Hence one must conclude it isn't internal rounding in the sin() function itself causing the difference in the original code..

As before, I don't have time to really dig, but can show some results that are expected to compare against...

Answer 3

[Djamil Boulahbal]

Thank you to all who contributed. I truly appreciate that. The chirp function did the trick for me.

I'm now rather curious at the inner workings of the chirp function, ... but that's for another day.

Cheers