Creating a sinewave with logarithmic frequency

56 次查看（过去 30 天）
Djamil Boulahbal
Djamil Boulahbal 2024-10-24，20:48
评论： Djamil Boulahbal 2024-10-25，22:19
Ran in:
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
  4 个评论
显示 2更早的评论
dpb
dpb 2024-10-25，15:03
编辑：dpb 2024-10-25，17:48
Ran in:
% 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
Djamil Boulahbal 2024-10-25，19:49
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.

请先登录，再进行评论。

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

采纳的回答

Star Strider
Star Strider 2024-10-24，23:42
Ran in:
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.
.
  5 个评论
显示 3更早的评论
Star Strider
Star Strider 2024-10-25，21:10
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
Djamil Boulahbal 2024-10-25，22:19
Thanks again, you've been very thorough and helpful.

请先登录，再进行评论。

更多回答（2 个）

Djamil Boulahbal
Djamil Boulahbal 2024-10-25，0:19
编辑：Djamil Boulahbal 2024-10-25，0:27
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?
  1 个评论
dpb
dpb 2024-10-25，16:56
Ran in:
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...

请先登录，再进行评论。

Djamil Boulahbal
Djamil Boulahbal 2024-10-25，19:54
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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by