Creating a sinewave with logarithmic frequency

2024 10 24
3 个回答
回答已采纳
更新时间：2024 10 25
17 次查看（30 天）
Ran in:

0 个投票

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更早的评论 隐藏 2更早的评论

dpb
dpb 2024-10-24
Ran in:
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;
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
Djamil Boulahbal 2024-10-25
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
dpb 2024-10-25
编辑：dpb 2024-10-25
Ran in:
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
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
Ran in:

0 个投票

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更早的评论 隐藏 3更早的评论

Djamil Boulahbal
Djamil Boulahbal 2024-10-25
编辑：Djamil Boulahbal 2024-10-25
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
Star Strider 2024-10-25
Ran in:
Ran in:
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
Djamil Boulahbal 2024-10-25
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
Star Strider 2024-10-25
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
Thanks again, you've been very thorough and helpful.

请先登录，再进行评论。

更多回答（2 个）

Djamil Boulahbal
Djamil Boulahbal 2024-10-25
编辑：Djamil Boulahbal 2024-10-25

0 个投票

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 个评论
显示 -1更早的评论 隐藏 -1更早的评论

dpb
dpb 2024-10-25
Ran in:
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

0 个投票

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

类别

在 帮助中心 和 File Exchange 中查找有关 Title 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by