21 views (last 30 days)

Show older comments

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.

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.

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

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

Start Hunting!