Sinusoidal signal with variable instantaneous frequency

11 次查看（过去 30 天）
Sowmya MR
Sowmya MR 2017-6-12
编辑： John BG 2017-6-21
Below is a simple matlab code that generates a sinusoidal signal with user specified no of cycles. How do i change this such that each cycle has different instantaneous frequency? For example, if i have 10 cycles, then there should be 10 random frequencies in total where each cycle has different frequency.
clc
clear
f1=120;
fs=f1/60; %sampling frequency depends on f1
no_of_cycles=10;
f=no_of_cycles/60; %Frequency of sinusoid
duration=1;%%duration of the sinusoidal signal in minutes
A=1;
t=1/fs:1/fs:duration*no_of_cycles*1/f; %time index
y1=A*sin(2*pi*f*t);%%simulated sinusoidal signal

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

采纳的回答

John BG
John BG 2017-6-12
编辑：John BG 2017-6-12
Hi Sowmia MR
This is John BG ( <mailto:jgb2012@sky.com jgb2012@sky.com> )
1.
One way to answer your question is, without losing constant amplitude, which is a good quality of FM:
f1_range=120*[1:1/3:4]
f1=f1_range(randi([1 9],1,9))
A=1;
y1=0
for k=1:1:numel(f1)
y1c=A*sin(2*pi/f1(k)*[1:1:f1(k)]);
y1=[y1 y1c];
end
plot(y1)
.
2.
Another way would be constraining the amount of samples for each frequency to a constant value, that would be of 1 cycle of f0=120Hz
f0=120
f1_range=f0*[.7:.1:1.3]
f1=f1_range(randi([1 numel(f1_range)],1,numel(f1_range)))
A=1;
y1=0
for k=1:1:numel(f1)
y1c=A*sin(2*pi*f1(k)/f0*[0:1/f0:1]);
y1=[y1 y1c];
end
plot(y1)
.
3.
Note that now it's not 10 cycles but 7, 7*120 are the 840 samples shown above 2nd graph.
Recommended reading 'Telecommunication Systems' by Bruce Carlson
<http://een.iust.ac.ir/profs/Beheshti/communication%20system/Carlson%20-%20Communication%20Systems%20-%204ed%20.pdf>
4.
Please not that f1_range is on purpose chosen 'not too wide'.
If so alias occurs, data loss, the output signal and plot may be meaningless.
5.
If you want to learn more about FM basics with MATLAB, please check my accepted answer to question in this link:
https://uk.mathworks.com/matlabcentral/answers/318969-fm-modulation-without-using-fmmod
there's an explanation how to FM modulate a signal, with MATLAB and with SIMULINK
So, Sowmia MR
if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?
To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link
thanks in advance
John BG
<mailto:jgb2012@sky.com jgb2012@sky.com>
  4 个评论
显示 2更早的评论
John BG
John BG 2017-6-12
编辑：John BG 2017-6-13
Yes, it can be done, you want 120 samples in each cycle, AND amplitude continuity.
One can achieve such amplitude continuity, with a phase carry:
clear all;clc;clf
f0=120
f1_range=120*[1:1/3:4]
f1=f1_range(randi([1 9],1,9))
A=1;
y1=0
dph=0
for k=1:1:numel(f1)
y1c=A*sin(2*pi*f1(k)/f0*[0:1/f0:f1(k)/f0-1/f0]+dph);
y1=[y1 y1c([1:1:f0])];
dph=floor(f0*rem(f1(k),f0)/f0)+dph
end
plot(y1)
the random sequence of frequencies of this example is
f1 =
240.0000 320.0000 160.0000 360.0000 120.0000 200.0000 120.0000 120.0000 400.0000
.
Repeat to test different f1 patterns.
Note that sometimes there is a small gl
itch
It's caused by the rounding of the command floor and the simplified way the carry is implemented.
You may find also of interest the Hull University introduction to GMSK, the modulation used in GSM, document available here:
http://www.emc.york.ac.uk/reports/linkpcp/appD.pdf
GSM uses the modulation GMSK that departs from this basic FM (or FSK) precisely to avoid constellation zero crossings (not to be confused with base band signal zero crossings or modulated signal zero crossings) avoiding sharp transitions of the VCO input signal using a gaussian filter.
GMSK also uses In phase and Quadrature, that is more robust than just using one real carrier.
Regards
John BG
<mailto:jgb2012@sky.com jgb2012@sky.com>
John BG
John BG 2017-6-16
编辑：John BG 2017-6-21
removing command floor
clear all;clc;clf
f0=120
f1_range=120*[1:1/3:4]
f1=f1_range(randi([1 9],1,9))
A=1;
y1=0
dph=0
for k=1:1:numel(f1)
y1c=A*sin(2*pi*f1(k)/f0*[0:1/f0:f1(k)/f0-1/f0]+dph);
y1=[y1 y1c([1:1:f0])];
dph=rem(f1(k),f0)+dph
end
plot(y1)
fixes the phase glitches shown in the last graph of this answer: no floor better.
John BG
<mailto:jgb2012@sky.com jgb2012@sky.com>

请先登录，再进行评论。

更多回答（1 个）

Walter Roberson
Walter Roberson 2017-6-12
Figure out how many samples there are per "cycle" for the purpose of your "10 cycles". The last time for each cycle should end 1 sample before the next full cycle boundary. So for example 10 samples per "cycle" should correspond to relative times 0, 1/10, 2/10, 3/10, 4/10, 5/10, 6/10, 7/10, 8/10, 9/10, and should not go all the way to 1 (which would be 11 samples). If you end up using linspace along the way to generate N points, then use linspace(first_time, last_time, N+1 ) and then throw away the last sample.
Now, generate the list of instantaneous frequencies for each cycle. You should probably limit those to being no more than 1/2 the sampling frequency, to avoid aliasing.
instant_frequency_for_cycle = fs/2 * rand(1, no_of_cycles)
Now replicate each of those by the number of samples per cycle:
f = reshape( repmat( instant_frequency_for_cycle(:), 1, samples_per_cycle ), 1, [])
or if you have a sufficiently new MATLAB,
f = repelem( instant_frequency_for_cycle, samples_per_cycle )
or in any version you can use
f = kron(instant_frequency_for_cycle, ones(1, samples_per_cycle) )
After that it is a simple modification to what you had:
y1 = A * sin(2*pi*f.*t); %%simulated sinusoidal signal
the f*t got changed to f.*t . With f and t being the same length, this gives an instant frequency to every sample.
  3 个评论
显示 1更早的评论
Sowmya MR
Sowmya MR 2017-6-12
Hi John, Thanks for the great insight. Could you please help me with changing my code to work for any cycles?
Walter Roberson
Walter Roberson 2017-6-12
John, you missed my "With f and t being the same length". My outline does not care about the 120 of the original code. I did not assume that "cycle" for the purpose of "10 cycles" corresponds to "hertz".

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Transforms 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by