Frequency sweeping in a sinusoidal signal

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 

Answer 1

[William Rose]

@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.

[William Rose]

I just noticed that the code above does not initialize p, x, and y to their final size. The reason is that, on the very last pass, when i=j, the code assigns values to the element j+1 of each vector. Therefore the vectors need to get bigger by 1, on that pass. On all previous passes, they are already big enough, due to their initialization. Therfore it would be better to change the lines in the code above to

p=zeros(1,j+1);
x=zeros(1,j+1);
y=zeros(1,j+1);

so that the vectors do not need to be enlarged on the final loop pass.

[Nneka Onubogu]

@William Rose

Thank you again for your reply and thanks for teaching me how to initialize vectors before running a loop.

When I use dt =4e-6 in the first section of my code, after I run it, I get ‘bb’ as a complex double so I am not able to plot a graph.

In the second section of my code, I tried your new suggestion. The code runs and stops at this error message:

In an assignment A(:) = B, the number of elements in A and B must be the same.

Error in pumpmodunew2 (line 58)

p(i)=(2.7324*10^-22)*pinitial*cp*i;

I think its because in the example you used here, you assumed cp = 1, which is a scalar. But cp is actually a vector:

cp=A0+A1*sin(theta); %signal vector

y1(i) is the white noise defined in section 1 of the code.

[William Rose]

@Nneka Onubogu, I don;t know why bb is complex. Please post the current version of your code so I can see all of it and run it.

If you change line 58 to the following,

p(i)=(2.7324*10^-22)*pinitial*cp(i)*i;

you should no get an error. The factor i on the right hand side will act as a ramp function, as the loop goes from i=1 to imax. Therefore, when the loop has completed, p will be a chirp function whose amplitude increases linearly with time.

I recommend changing the variable name j to imax. j is often used as the index variable for a loop, rather than the termination value for a loop over i. Using imax will make the code more readable to others and will avoid confusion with .

An alternative way to compute vector p is to do it outside the loop, before the loop:

p=2.7324E-22*pinitial*cp.*(1:imax)

The line above multiplies the vector cp by the vector (1:), using the dot-times operator, ".*", which multiplies two vectors or arrays element-by-element, giving a vector or array with the same size as the inputs. It will be an error if the array sizes are unequal. If you use the line above, then you do not need to initialize p to a vector of zeros, because it will not be growing inside a loop.

Change the loop maximum value from imax to imax-1:

for i=1:imax-1
    %x(i+1)=...; 
    %y(i+1)=...;
end

so that the vectors x and y have length imax, rather than imax+1, when the loop is done.

The lines in the loop indicate that you are calculating x and y by integrating two differential equations:

I recommend that you confirm that those are the correct differential equations.

[Nneka Onubogu]

Thanks for your in-depth explanation. This is the current version of the code (section 1 and section 2). I have changed line 58 accordingly and j to imax. I still get an error message with the changes made.

 clc; clearvars; warning off;

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;
%dt=4e-6;  %time step (s), should be at least 5x smaller than 1/fc2
imax=1*10^8;%iterations for accuracy
%j=2500000;
fs=1/dt; %sampling frequency in Hz
t=(0:imax)*dt;     %time vector
%t=(0:j)*dt/2.4156e8;
p=zeros(1,imax+1);
x=zeros(1,imax+1);
y=zeros(1,imax+1);
y1 = wgn(imax,1,-35,50); %white noise 
for k=20.5 % pumping power
    pinitial=(k/10)*1e19; %pump power 
end    
for i=1:imax
   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 


%%

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=0.4;       %amplitude of oscillation
Tend=30;     %chirp duration (s)
t=0:dt:Tend;                                %time vector
theta=2*pi*(fc1*t+(fc2-fc1)*t.^2/(2*Tend)); %theta vector
cp=A0+A1*sin(theta);                         %signal vector
figure;
subplot(3,1,1); plot(t,cp,'r');              %plot entire signal
subplot(3,1,2); plot(t(1:1000),cp(1:1000),'r');  %plot beginning part
subplot(3,1,3); plot(t(end-100:end),cp(end-100:end),'r'); %plot end part
figure;
spectrogram(cp,1024,512,1024,1/dt,'yaxis');  %plot spectrogram

u=0; % amplitude (2nd order)
f2amp=u/100;    
r=40; % amplitude (1st order)
f1amp=r/100;
%p=2.7324E-22*pinitial*cp.*(1:imax);
   for i=1:imax
       p(i)=(2.7324*10^-22)*pinitial*cp(i)*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 

The correct differential equations are:

Thanks a lot for all your help.

[William Rose]

@Nneka Onubogu,

In the second part of the code, you compute x(t) and y(t). The differential equation for y(t) involves p(t). Signal p(t) is the chirp signal modulated by a ramp in amplitude. I am getting an error in this part because the duraitons of x(t) and y(t) do not match the duration of cp(t) and p(t). You have set the chirp duration to be 30 seconds, by choosing Tend=30 in the chirp code. But the duraion of x(t) and y(t) is the product of dt=1e-9 and imax=1e8. Therefore the duration of x(t) and y(t) is 0.1 seconds. Which is correct? They must match.

You generate white Gaussian noise in section 1 of your code. The power of this noise is inversely proportional to the step size, so there is ten times more power when dt=1e-9 than when dt=1e-8. Threfore the step size choice could influence the simulaiton results significantly. This is not good, because real physical systems do not have a discrete step size. Therefore I encourage you to considet whether the physical noise which you are simulating with gwn() has an upper frequency limit. Then we can generate band-limited Gaussian white noise. The step size shoudl be 4-5 times smaller than the period of the highest frequency noise which you wish to simulate. In your code as it stands, there is no upper frequency limit to the noise.

When simulating chaotic dynamics, you do not need to add external noise. SInce the chaotic behavior can look like noise, many investigators will first simulate with no noise, to understand the chaotic dynamics without noise, and then add noise. You may wish to condsider this option.

I obtained your paper, "Dynamics of a highly sensitive Erbium-doped Fabry-Perot Fiber (EDFBF) Laser Sensor under pump modulation", from the 2019 Microoptics Conference in Japan. It is very nice. I wish there was a time scale shown on the x-axis of the traces in figure 3, in order to uderstand the likely duration for the chirp and the simulation.

[William Rose]

Another way of writing the differential equations:

This is easier for me to read. This helps me see that the rate constants in the diff.eq. for x are . The values are approximately 9.2e+6, 2.6e-6, 2.0e-6, and 3.8e+6.. The units are inverse seconds, if x and y are dimensionless. The rate constants vary by 12 orders of magnitude. It means that some system processes occur at the microsecond time scale and some parts occur at the million-second time scale. I hope this is an error. If it is correct, then the system is very, very stiff, and you will need to simulate millions of seconds, with sub-microsecond resolution, to understand its behavior.

The rate constants in the diff.eq. for y are . The values are approximately 3.8e+6, 1.3e+0, 1.9e+6. The units are inverse seconds, if x and y are dimensionless. This is not as bad as the differential equations for x, but it is still a very stiff system, with rate constants ranging over six orders of magnitude. The simulation time step size should be smaller than the reciprocal of the fastest rate constant.

Check the units of variables and constants on both sides of differential equations for consistency. What are the units for the constant -1 on the right hand side of the diff.eq. for dy/dt? If -1 is dimensionless, then it must be the case that y has units of seconds, because dy/dt has units of [y units]/seconds. Then what are the units for x?

[William Rose]

In your code, you trned off the warning messages at the top. I woudn;t do this. Warning messages can be useful. I ran your code, all the way through, with dt=1e-9 seconds, and with chirp duration=0.1 seconds, to match the duration of x(t) and y(t). Even with dt=1e-9, I get a complex array for bbb. I know this because I disabled the suppression of warnings, and I see this warning message:

>> chirpGenerator2
x(t): mean=0.000,sd=0.000, min=0.000, max=0.002
y(t): mean=2.373,sd=0.206, min=0.000, max=2.481
Part 1 done.
Chirp done, now start part 2.
Warning: Imaginary parts of complex X and/or Y arguments ignored 
> In chirpGenerator2 (line 105) 
Part 2 done.

I have formatted my console display as code. This causes some strange colors, but is otherwise OK. The code I ran is attached.

[William Rose]

I suspect we are getting complex values for bbb because the vector xx, whose spectrum we are computing, has NaNs in it. This happens because x has NaNs in it. This happens because the differential equations for x and y in part 2 lead to crazy results. I obtained similar reults with your version of the differential equations and with my simplified version of the differential equations. The results are not exactly the same because the random noise is different on each run. Here is my current version of the code. It displays staistics for x(t) and y(t) after part 1 and after part 2. It takes about 2 minutes to run on my computer. It produces the console output below.

>> chirpGenerator201
x(t): mean=0.000,sd=0.000, min=0.000, max=0.002
y(t): mean=2.373,sd=0.206, min=0.000, max=2.481
Part 1 done.
Chirp done, now start part 2.
x(t): mean=NaN,sd=NaN, min=-Inf, max=262.678
y(t): mean=NaN,sd=NaN, min=-1044.955, max=Inf
Warning: Imaginary parts of complex X and/or Y arguments ignored 
> In chirpGenerator201 (line 113) 
Part 2 done.

[William Rose]

In part 1, x(t) and y(t) have reasonable values, and the power spectrum is real. In part 2, x(t) and y(t) have off-scale (NaN) values, and the power spectrum is complex. The only difference between parts 1 and 2 is that p(t) is defined differently. In part 1, p(t) is a constant. In part 2, p(t) is a chirp signal multiplied by an amplitude ramp.

I changed the definition of p(t) in part 2 by removing the amplitude ramp component. This fixes the problem with NaNs: now there are no NaNs in x(t) or y(t). It also fixes the problem of a complex power spectrum: now the power spectrum is real. The chirp in part 2 goes from 1 kHz to 30 kHz in 0.1 seconds.

Here are plots from the code. Notice the difference in the power spectrum from 1 to 30 kHz in part 1 and part 2.

The console outout is below:

>> chirpGenerator203
p(t): mean=0.001,sd=0.000, min=0.001, max=0.001
x(t): mean=0.000,sd=0.000, min=0.000, max=0.002
y(t): mean=2.373,sd=0.205, min=0.000, max=2.480
Part 1 done.
Chirp done, now start part 2.
p(t): mean=0.000,sd=0.002, min=-0.002, max=0.002
x(t): mean=0.000,sd=0.000, min=0.000, max=0.000
y(t): mean=0.685,sd=0.089, min=0.000, max=1.283
Part 2 done.

The code is attached.

Minor question: The constant that multiplies pinitial, when calculating p, is different by a factor of 10 in part 1 (2.73e-23) and part 2 (2.73e-22). Is this correct?

[Nneka Onubogu]

@William Rose

Thanks a lot for all your comments.

If i correctly understand what you mean by ramp, I think signal p(t) Is modulated by a ramp in frequency not amplitude, the amplitude is constant. I initially set the duration of x(t) and y(t) as 0.1 seconds and then changed it to 30 secs after I consulted you. I think I prefer making it 30 secs, so that it is 1 sec per frequency.

About the Gaussian White noise, please can you give me an example of step size? Are you referring to this: Add white Gaussian noise to input signal - MATLAB (mathworks.com)? I will check whether the physical noise which I am simulating with gwn() has an upper frequency limit; I am not sure it has. I previously simulated without white noise and the results were not as expected, so I thought it was proper to include white noise as there was little background noise during the experiments. I read online that band-limited white Gaussian noise cannot be generated but band-limited Gaussian noise as white noise means the power spectral is flat. To add a band-limit to my present Gaussian white noise, can I do it like this: How can I generate band-limited Gaussian white noise?? - MATLAB Answers - MATLAB Central (mathworks.com)? As you have suggested, I will not include white noise first and after analysing my results, I will see if there is need to include it. Thanks for your comment on my conference paper. I included a time-scale (ms) in my other related journal papers on Optik.

The way you have written the equations Is much better for clarity. The x and y are dimensionless (x and y are the normalized laser power density and inversion population of my laser model). Rate equations are often a so-called “stiff” system of differential equations, involving very different time constants. The constant -1 on the right hand side of the diff.eq. for dy/dt is dimensionless. The constant that multiplies pinitial, when calculating p, is actually (2.73e-23), I had changed it to (2.73e-22) to see what results I get.

Thanks a lot for all you efforts and thanks for the plots and changes made in my code. I will go through it and comment on the results.

[William Rose]

@Nneka Onubogu, will follow up by email.

[Nneka Onubogu]

@William Rose, ok. Thank you very much. I will accept your answer as you have actually shown me how to generate the chirp frequency.

Answer 2

[William Rose]

@Nneka Onubogu,

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]

@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.

Answer 3

[William Rose]

@Nneka Onubogu,

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 forabove. 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.

Answer 4

[Nneka Onubogu]

@William Rose

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.