# Frequency sweeping in a sinusoidal signal

21 views (last 30 days)
Nneka Onubogu on 28 May 2021
Commented: Nneka Onubogu on 7 Jun 2021
I want to generate a chirp signal by sweeping my frequency from 0 kHz to 30 kHz. Can i achieve it this way?
a1=2.4;
a2=6.9*10^-13;
a3=5.1*10^-13;
b1=3.5*10^-7;
b2=2.6*10^-7;
b3=0.5;
dt=10^-9;
x(1)=0;
y(1)=0;
j =1*10^8;%iterations for accuracy
fs=1/dt; %sampling frequency
t=(0:j)*dt/2.4156e8;
y1 = wgn(j,1,-35,50); %white noise
for k=20.5 % pumping power
for i=1:j
pinitial=(k/10)*1e19; %pump power
p(i)=(2.7324*10^-23)*pinitial*(1);
x(i+1)=(((x(i)*y(i))/b2)-(((a1+y1(i))/b2)*x(i))+((a2/b2)*y(i))+(a3/b2))*dt+x(i);
y(i+1)=(-((x(i)*y(i))/b2)-((b1/b2)*y(i))-1+((p(i)/b2)*(1-b3*(1))))*dt+y(i);
end
[bb,aa]=pwelch(x./9.8317e-17,[],[],[],fs); % aa is frequency, bb is amplitude
figure
semilogy(aa(1:3357),bb(1:3357))
grid minor
end
for fp = 1000+((30000-1000)/(j*dt/3.33e-8))*t(1:end-1)
for u=0 % amplitude (2nd order)
f2amp=u/100;
for r=40 % amplitude (1st order)
f1amp=r/100;
for i=1:j
p(i)=(2.7324*10^-22)*pinitial*(1+(f1amp*sin(2*pi*dt*i.*fp/1))+(f2amp*sin(2*pi*dt*i.*fp/2)));
x(i+1)=(((x(i)*y(i))/b2)-(((a1+y1(i))/b2)*x(i))+((a2/b2)*y(i))+(a3/b2))*dt+x(i);
y(i+1)=(-((x(i)*y(i))/b2)-((b1/b2)*y(i))-1+((p(i)/b2)*(1-b3*(1))))*dt+y(i);
end
end
end
end
downsampling=1;
xx=downsample(x./9.8317e-17,downsampling);
[bbb,aaa]=pwelch(xx,[],[],[],fs/downsampling);
figure()
subplot(2,1,1)
plot(t,xx) % time domain
grid minor
subplot(2,1,2)
semilogy(aaa(1:3357),bbb(1:3357))% frequency domain
grid minor

William Rose on 31 May 2021
@Nneka Onubogu, You can initialize the vectors p, x, and y before the loop as follows:
dt=4e-6; %(s)
fs=1/dt; %(Hz)
j=25000; %make this 250K if you want 1 sec duration
t=(0:j)*dt; %time vector
cp=1; %Define cp before the loop.
p=zeros(1,j); %or p=5*ones(1,j) if you want a vector of all fives, etc.
x=zeros(1,j);
y=zeros(1,j);
pinitial=2.05e19;
a1=1; a2=1; a3=1; b1=1; b2=1; b3=1; %define these constants as desired
for i=1:j
p(i)=(2.7324*10^-22)*pinitial*cp*i;
x(i+1)=(((x(i)*y(i))/b2)-(((a1+y(i))/b2)*x(i))+((a2/b2)*y(i))+(a3/b2))*dt+x(i);
y(i+1)=(-((x(i)*y(i))/b2)-((b1/b2)*y(i))-1+((p(i)/b2)*(1-b3*(1))))*dt+y(i);
end
%downsampling=1;
%xx=downsample(x./9.8317e-17,downsampling);
%Downsampling by 1 does nothing. When you divide a vector by a constant,
%you do not need to "dot-divide". Regular division is fine.
xx=x/9.8317e-17;
[bbb,aaa]=pwelch(xx,[],[],[],fs);
figure
subplot(2,1,1)
plot(t,xx) % time domain
grid minor
subplot(2,1,2)
semilogy(aaa,bbb)% frequency domain
grid minor Code above runs without error. I changed y1(i) to y(i) in the line x(i+1)=... because y1(i) was probably a mistake.
Nneka Onubogu on 7 Jun 2021
@William Rose, ok. Thank you very much. I will accept your answer as you have actually shown me how to generate the chirp frequency.

William Rose on 28 May 2021
If you highlight your code in your posting, and then click the "code" icon at the top of the pane, it will format your code nicely, and more importantly, it will allow others to run your code. Therefore please do this - and check that your code is in fact runnable, or if it isn't, tell us what error you get.
For generating a chirp and other signals with varying frequency, the built-in vco() function (here) is very convenient.
Good luck.
Nneka Onubogu on 28 May 2021
@William Rose thank you so much for your comment and thanks for teaching me on how to format my code as I am new to Matlab. I will do that right now and will post the complete code so that it can be run by others.
When I run the second section of my code, it takes forever and after i get a result, my (t,xx) plot is not a chirp signal but rather a signal for only one frequency.

William Rose on 29 May 2021
Thak you for reformatting the code part of your original post.
The original post asks how do I do a chirp from 0 to 30 kHz, but there is a lot more going on here than just a chirp.
Matlab has a chirp() command and a vco() command. vco() is more general than chirp, but can be used to make a chirp.
Your code will crash my machine if I try to run it, or it will take forever, because you have a for loop that is set to run 100 million times. Each time, the x and y arrays grow by one element. By the end, x() and y() will have 100 million elements each. That will never work. You shoud initialize arrays before you run the loop, because otherwise, each array gets copied to a new array on each loop pass, which gets very slow when the array is huge. And arrays that huge will probably not work anyway.
You also have a lot of weird constants. For example:
t=(0:j)*dt/2.4156e8;
Why do you divide by 241 million? If you define the t vector this way, then it will not have the spacing "dt" that you have specified previously.
Another example of weird constants:
pinitial=(k/10)*1e19; %pump power (k=20.5 defined previously)
p(i)=(2.7324*10^-23)*pinitial*(1);
In what possible units is the pump power 2x10^19? In the next line, you multiply that gigantic number by a very tiny number, resulting in p(i)=0.0005 (approximately). Why use such extreme and offsetting constants? Also, p(i) defined this way in this loop is the same on every loop iteration, so you should define it once, outside the loop, to save time.
Other things:
for k=20.5, % pumping power ...
end
for u=0, % amplitude (2nd order) ...
end
for r=40, % amplitude (1st order) ...
end
The for loops above have only a single value for the loop variable, so each loop will only execute once. This will run, but it does not make sense to set up for loops that only run once.
I suspect the following code is where you try to generate a chirp. I have added a few comments.
for fp = 1000+((30000-1000)/(j*dt/3.33e-8))*t(1:end-1)
for u=0 % amplitude (2nd order) (This will only run once)
f2amp=u/100; %(f2amp will equal zero, so why bother?)
for r=40 % amplitude (1st order) (this will only run once)
f1amp=r/100; %(f1amp=0.4, so define it outside the loop)
for i=1:j
p(i)=(2.7324*10^-22)*pinitial*(1+(f1amp*sin(2*pi*dt*i.*fp/1))+(f2amp*sin(2*pi*dt*i.*fp/2)));
x(i+1)=(((x(i)*y(i))/b2)-(((a1+y1(i))/b2)*x(i))+((a2/b2)*y(i))+(a3/b2))*dt+x(i);
y(i+1)=(-((x(i)*y(i))/b2)-((b1/b2)*y(i))-1+((p(i)/b2)*(1-b3*(1))))*dt+y(i);
end
end
end
end
It appears that p(i) is the variable that is supposed to chirp, and the chirp frequency is supposed to vary from 1000 to 30000. The code above will not work, for various reasons which would take a while to explain.
Let's go back to the beginning. The min and max chirp frequencies are and . You specify dt=1e-9, but you can use a much longer step, which will allow many fewer samples and therefore faster execution and less risk of memory overflow. The fastest frequency is 30 kHz, which has a perod of 1/(30kHz)=33 microseconds. Therefore, if we choose dt=4 microseconds, we will have about 8 samples per cycle at the highest frequency, which is enough, and more samples than that at the lower frequencies of the chirp. So let's specify dt=4e-6.
We want frequency f to equal fc1 at t=0 and f=fc2 at t=Tend. We can write the frequency function as follows: How long should the chirp take? Let's try one second, i.e. . This is not a crazy duration, since it only takes 1 msec to get through 1 cycle at the slowest frequency.
To make a chirp, you need a signal such as where the rate of change of is steadily increasing. ( is the mean value of p and is the amplitude of the sinusoidal part.) For a regular un-chirped sine wave, the rate of change of is constant: , from which it follows that .
For a chirp, f is not constant. It is given by the equation for above. Therefore we have This is a simple and solvable differential equation: Integrate both sides, and assume when t=0, and you get Therefore the Matlab code to make a chirp from 1 kHz to 30 kHz, lasting one second, is
fc1=1000; %start frequency (Hz)
fc2=30000; %end frequency (Hz)
dt=4e-6; %time step (s), should be at least 5x smaller than 1/fc2
A0=0; %mean value of signal
A1=1; %amplitude of oscillation
Tend=1; %chirp duration (s)
t=0:dt:Tend; %time vector
theta=2*pi*(fc1*t+(fc2-fc1)*t.^2/(2*Tend)); %theta vector
p=A0+A1*sin(theta); %signal vector
figure;
subplot(3,1,1); plot(t,p,'r'); %plot entire signal
subplot(3,1,2); plot(t(1:1000),p(1:1000),'r'); %plot beginning part
subplot(3,1,3); plot(t(end-100:end),p(end-100:end),'r'); %plot end part
figure;
spectrogram(p,1024,512,1024,1/dt,'yaxis'); %plot spectrogram
The code generates 2 figures. The first figure shows the whole signal, the first 4 milliseconds, and the last 0.4 milliseconds. This figure show that frequency is 1 kHz at the start and 30 kHz at the end, as desired. The second figure is the spectrogram, or time-dependent power spectrum. It shows that the frequency increases from 1 kHz to 30 kHz from t=0 to t=1 s.  Change Tend, A0, or A1 to alter the chirp duration, mean, or amplitude.

Nneka Onubogu on 31 May 2021
Thank you very much for taking your time to look into my code and making corrections line by line. I am very grateful. Actually, those weird constants were gotten from some equations based on my laser set-up and rate equations for laser. The pump power of 2*10^19 gives a resonance frequency of 10 kHz (when you plot the graph, you can see a peak at 10 kHz), so I change the value of ‘k’ to be able to get different resonance frequencies. I forgot to indicate that the code is in two parts. I actually run the first part of the code that stops at the “end” before the ‘fp’ for loop.
Please how can I initialize arrays before running the loop? Please can you show me an example?
I will remove unnecessary for loops. Thanks for the advice and corrections. Thanks also for explaining in details on how to obtain the chirp and for plotting it out to show me too.
The main aim of my simulation is to obtain the dynamic behavior of my fiber laser system. To do this, I needed to sweep my modulation frequency (0 to 30kHz) so as to get the time domain which should be the chirp signal and then measure the maximum peak amplitude at each frequency from which I will plot a bifurcation diagram. That Is the essence of this part of my code:
for i=1:j
p(i)=(2.7324*10^-22)*pinitial*(1+(f1amp*sin(2*pi*dt*i.*fp/1))+(f2amp*sin(2*pi*dt*i.*fp/2)));
x(i+1)=(((x(i)*y(i))/b2)-(((a1+y1(i))/b2)*x(i))+((a2/b2)*y(i))+(a3/b2))*dt+x(i);
y(i+1)=(-((x(i)*y(i))/b2)-((b1/b2)*y(i))-1+((p(i)/b2)*(1-b3*(1))))*dt+y(i);
end
end
end
end
downsampling=1;
xx=downsample(x./9.8317e-17,downsampling);
[bbb,aaa]=pwelch(xx,[],[],[],fs/downsampling);
figure()
subplot(2,1,1)
plot(t,xx) % time domain
grid minor
subplot(2,1,2)
semilogy(aaa(1:3357),bbb(1:3357))% frequency domain
grid minor
Based on your corrections, would it be correct to change the above section of my code in this way:
Assuming I call the signal vector cp,
for i=1:j
p(i)=(2.7324*10^-22)*pinitial*cp*i;
x(i+1)=(((x(i)*y(i))/b2)-(((a1+y1(i))/b2)*x(i))+((a2/b2)*y(i))+(a3/b2))*dt+x(i);
y(i+1)=(-((x(i)*y(i))/b2)-((b1/b2)*y(i))-1+((p(i)/b2)*(1-b3*(1))))*dt+y(i);
end
downsampling=1;
xx=downsample(x./9.8317e-17,downsampling);
[bbb,aaa]=pwelch(xx,[],[],[],fs/downsampling);
figure()
subplot(2,1,1)
plot(t,xx) % time domain
grid minor
subplot(2,1,2)
semilogy(aaa(1:3357),bbb(1:3357))% frequency domain
grid minor
When I run it, I get an error on the p(i) line. I want to be ale to plot :
plot(t,xx) % time domain and semilogy(aaa(1:3357),bbb(1:3357))% frequency domain.
Thanks a lot.