Question

When I do a spectrogram in matlab / octave I can create a swept signal that looks like the RED plot line below. But how can I create a swept signal like the BLUE line in the 1st plot using the equation below.

thanks to Daniel and David for getting me this far with the code is below

startfreq=200;
fs=44100;
endfreq=20;
dursec= 10;%duration of signal in seconds
t=(0:dursec*fs)/fs; %Time vector
alpha=log(startfreq/endfreq)/dursec;
sig = exp(-j*2*pi*startfreq/alpha*exp(-alpha*t));
sig=(sig/max(abs(sig))*.8); %normalize signal
wavwrite([sig'] ,fs,32,strcat('/tmp/del.wav')); %export file
specgram(sig,150,400);

1st plot Plot that I'm trying to reproduce


2nd plot Plot with 200hz-20hz works


How can I fix the the equation in the variable sig to get it to look like the BLUE line in the 1st plot?

3rd plot Plot with 20hz-200hz which doesn't do what I want


Was it helpful?

Solution

This question is almost an month old, so you might have figured this out by now. Here's an answer in case you are still interested.

It appears that your current model for the frequency is

freq(t) = b*exp(-alpha*t)

with

freq(0) = b = startfreq
freq(dursec) = b*exp(-alpha*dursec) = endfreq

There are two free parameters (b and alpha), and two equations. The first equation, b = startfreq, gives us b (trivially).

Solving the last equation for alpha gives

alpha = -log(endfreq/startfreq)/dursec
      = log(startfreq/endfreq)/dursec

So

freq(t) = startfreq * exp(-alpha*t)

To use this as the instantaneous frequency of a frequency-swept signal, we need the integral, which I'll call phase(t):

phase(t) = -(startfreq/alpha) * exp(-alpha*t)

The (complex) frequency-swept signal is then

sig(t) = exp(2*pi*j * phase(t))

The real part of this signal is

sig(t) = cos(2*pi*phase(t))

That explains your current code. To generate a chirp whose frequency varies like the blue curve, you need a different model for the frequency. A more general model than the one used above is

freq(t) = a + b*exp(-alpha*t)

The requirements at t=0 and t=dursec are

freq(0) = a + b = startfreq
freq(dursec) = a + b*exp(-alpha*dursec) = endfreq

That's two equation, but we now have three parameters: a, b, and alpha. I'll use the two equations to determine a and b, and leave alpha as a free parameter. Solving gives

b = (startfreq - endfreq)/(1 - exp(-alpha*dursec))
a = startfreq - b

Integrating the model gives

phase(t) = a*t - (b/alpha)*exp(-alpha*t)

alpha is an arbitrary parameter. Following the formula from the first model, I'll use:

alpha = abs(log(startfreq/endfreq))/dursec

The following is a complete script. Note that I also changed the use of exp(-j*2*pi*...) to cos(2*pi*...). The factor 0.8 is there to match your code.

startfreq = 20;
endfreq = 200;
fs = 44100;
dursec = 10;  % duration of signal in seconds
t = (0:dursec*fs)/fs;  % Time vector

if (startfreq == endfreq)
    phase = startfreq * t;
else
    alpha = abs(log(endfreq/startfreq))/dursec;
    b = (startfreq - endfreq)/(1 - exp(-alpha*dursec));
    a = startfreq - b;

    phase = a*t - (b/alpha)*exp(-alpha*t);
endif

sig = 0.8 * cos(2*pi*phase);

wavwrite([sig'] ,fs,32,strcat('del.wav'));  % export file
specgram(sig,150,400);
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top