Amplitude modulation in scilab

[M. Hamed]

Disclaimer: this may be a tiny bit off topic for this group.

Now DSP is one of the things that I'm planning to dive deep into at some point. So I know I'm going to get a bunch of those "Get a DSP book". But since I'm going to be working on an AM radio, I decided I would like to understand the mixing process better and demonstrate it in software in a burn before you learn way.

So I fired up Scilab and did the following:

1- create a time array
2- create a bunch of sinusoids at frequencies from 550khz o 1700 khz spacedby 20khz and random phase.
3- Sum all of these sinusoids to create a composite wave.
4- perform FFT on the composite wave.
5- Plot the results

The code is here, if anyone is interested.

https://docs.google.com/document/d/1EQp43CK_ooTF66am9pXCOu5k-vEK21wWJhsdSa6NvZQ/edit?usp=sharing

Now I have to tell you I was a bit surprised by the result. I started playing with both sampling rate and and found the results quite interesting.

This is a picture of the FFT plot with the sampling rate a bit above Nyquist at 4M samples/s (Max AM frequency is 1.7MHz)

https://docs.google.com/file/d/0B1hjUVk4ytmvdXBELTRwanBlcTA/edit?usp=sharing

Now you can see 58 sinusoids which is the correct number created in the time domain, but I expected them to have equal amplitudes. I assume the decreasing amplitude is an artifact of the low sampling rate.

When I increase the sampling rate 10 times more and decrease the time to account for memory shortage, I got what I wanted. 58 sin waves with equal amplitude.

https://docs.google.com/file/d/0B1hjUVk4ytmvMmdvb3JwOGpneDQ/edit?usp=sharing

So far my signal sampling time window was 1 ms. As I decrease this further,I start getting more noise in the amplitudes of the sin waves. For examplehere with a time of 100 us.

https://docs.google.com/file/d/0B1hjUVk4ytmvWGRTdzBfRllZeTQ/edit?usp=sharing

It eventually reaches a point where the sinusoids are completely messed up,even at only half the previous time window 50 us and even with a 400M samples per second rate, much higher than nyquist.

https://docs.google.com/file/d/0B1hjUVk4ytmvWFlvWHZqQ2VJWms/edit?usp=sharing

I am sure there is a DSP explanation for all this. It gives me some motivation to study further. One possibility with relatively short time windows isthat you don't have enough time to capture all the information in the signal that it can't yet be approximated with a sin wave in the discreet domain.. Another possibility is that may be there has to be some relation that needs to be maintained between time window and sampling rate. Don't know.

Next on my plate is to do actual mixing with a sinusoids then do some filtering then another step of down-conversion and demodulation. I have no idea how to do digital filtering in scilab yet but I'll find out.

It was interesting to me, thought It's worth sharing and may be I can get some hints.

 

[Fred Abse]

The code is here, if anyone is interested.

<URL snipped>

Why not post it to this newsgroup, instead of that Javascript ridden hell?

 

[M. Hamed]

Thanks for the tips on windowing. I'll take look into it.

Why not post it to this newsgroup, instead of that Javascript ridden hell?

You are absolutely right. I will look into some other alternative that shares files as files.

Here is the code:

///////////////////////
// Program Constants

start_freq = 550e3; // starting frquency = 550 kHz
end_freq = 1700e3; // end frequency = 1700 kHz
freq_spacing = 20e3; // channel spacing = 20 kHz
t_end = 1e-2; // last time sample = 10 us
amplitude = 1;
sample_rate = 4e6;
random_phase = 1;
fftplot = 1;
nf = 20000 // number of frequency samples to plot

///////////////////////

// create time variable

time = 0:1/sample_rate:t_end;
szt = size(time, '*');

// create station frequencies spaced by 20khz
frq = start_freq:freq_spacing:end_freq;
// transpose
frq = frq';
szf = size(frq, '*');

//phase = (2*%pi/szf)2*%pi/szf):2*%pi;
phase = 2*%pi*rand(1,szf);
phase = phase'
phase_matrix = random_phase*phase*ones(1,szt);

// time frequency matrix, one row for each frequency and time is in columns
time_freq = amplitude*sin((2*%pi*frq*time)+phase_matrix);

// sum all the sinusoids
composite = sum(time_freq,;

// Now calculate and plot FFT
if fftplot then
N = szt;

y = fft(composite);

f=sample_rate*(0N/2))/N;
n=size(f,'*');

clf()
if nf=0 then
nf = n
end
plot(f(1:nf),abs(y(1:nf)))
end

//////////////////////////////////////////////////

 

[Fred Abse]

You are absolutely right. I will look into some other alternative that
shares files as files.

Here is the code:

Thanks.

Unfortunately it runs my scilab out of memory. I'll see whether I can
translate it into octave code.

 

[josephkk]

Disclaimer: this may be a tiny bit off topic for this group.

Now DSP is one of the things that I'm planning to dive deep into at somepoint. So I know I'm going to get a bunch of those "Get a DSP book". Butsince I'm going to be working on an AM radio, I decided I would like to understand the mixing process better and demonstrate it in software in a burn before you learn way.

So I fired up Scilab and did the following:

1- create a time array
2- create a bunch of sinusoids at frequencies from 550khz o 1700 khz spaced by 20khz and random phase.
3- Sum all of these sinusoids to create a composite wave.
4- perform FFT on the composite wave.
5- Plot the results

The code is here, if anyone is interested.

https://docs.google.com/document/d/1EQp43CK_ooTF66am9pXCOu5k-vEK21wWJhsdSa6NvZQ/edit?usp=sharing

Now I have to tell you I was a bit surprised by the result. I started playing with both sampling rate and and found the results quite interesting.

This is a picture of the FFT plot with the sampling rate a bit above Nyquist at 4M samples/s (Max AM frequency is 1.7MHz)

https://docs.google.com/file/d/0B1hjUVk4ytmvdXBELTRwanBlcTA/edit?usp=sharing

Now you can see 58 sinusoids which is the correct number created in the time domain, but I expected them to have equal amplitudes. I assume the decreasing amplitude is an artifact of the low sampling rate.

When I increase the sampling rate 10 times more and decrease the time toaccount for memory shortage, I got what I wanted. 58 sin waves with equal amplitude.

https://docs.google.com/file/d/0B1hjUVk4ytmvMmdvb3JwOGpneDQ/edit?usp=sharing

So far my signal sampling time window was 1 ms. As I decrease this further, I start getting more noise in the amplitudes of the sin waves. For example here with a time of 100 us.

https://docs.google.com/file/d/0B1hjUVk4ytmvWGRTdzBfRllZeTQ/edit?usp=sharing

It eventually reaches a point where the sinusoids are completely messed up, even at only half the previous time window 50 us and even with a 400M samples per second rate, much higher than nyquist.

https://docs.google.com/file/d/0B1hjUVk4ytmvWFlvWHZqQ2VJWms/edit?usp=sharing

I am sure there is a DSP explanation for all this. It gives me some motivation to study further. One possibility with relatively short time windows is that you don't have enough time to capture all the informationin the signal that it can't yet be approximated with a sin wave in the discreet domain. Another possibility is that may be there has to be some relation that needs to be maintained between time window and sampling rate. Don't know.

Next on my plate is to do actual mixing with a sinusoids then do some filtering then another step of down-conversion and demodulation. I have no idea how to do digital filtering in scilab yet but I'll find out.

It was interesting to me, thought It's worth sharing and may be I can get some hints.

Interesting. Much of what you are discovering is normal artifacts of the
FFT algorithms themselves. Or properties finite time FT. Yes, everything
that you saw in these examples. You could really profit from E. Oran
Brighams book on FFT and applications. Not much heavy math but much
explanation.

?-)

 

[josephkk]

<URL snipped>

Why not post it to this newsgroup, instead of that Javascript ridden hell?

Well, not everybody here has scilab. It does not help much for a google
grouper trying to post graphics.

Fred, if you want a wake up call try scilab.org itself with a recent
version of firefox.

Reccomendation for M. Hamed, get pan (news client and email client) and
use eternal-september for free posting. I think that they do carry
alt.binaries.schematic.electronic (abse or a.b.s.e in common parlance).

?-)

 

[Fred Abse]

Well, not everybody here has scilab. It does not help much for a google
grouper trying to post graphics.

I hardly ever use scilab, preferring octave, mainly because it makes a
better fist of running matlab stuff.
What I asked for wasn't graphics, but scilab code, ie. text.

Fred, if you want a wake up call try scilab.org itself with a recent
version of firefox.

No thanks, I climbed off the upgrade bandwagon aeons ago. Too many
dependency changes.

Reccomendation for M. Hamed, get pan (news client and email client) and
use eternal-september for free posting. I think that they do carry
alt.binaries.schematic.electronic (abse or a.b.s.e in common parlance).

Seconded.

 

[M. Hamed]

On Sat, 27 Apr 2013 16:55:56 -0700 (PDT), "M. Hamed"
Interesting. Much of what you are discovering is normal artifacts of the
FFT algorithms themselves. Or properties finite time FT. Yes, everything
that you saw in these examples. You could really profit from E. Oran
Brighams book on FFT and applications. Not much heavy math but much
explanation.

I'm currently trying to go through "The scientist and engineer's guide to dsp". Most of these DSP books try to include very little mathematics and I feel this takes away from my understanding. I am thinking I should go for a Signals and Systems book for a better groundwork then move to DSP (except that it takes a long time this way ).

Back in my school days we had a Signal Analysis course and the book was Signal Analysis by Papoulis. It was a bit difficult to follow back then but itwould probably be easier now. Can't find a good copy though.

I'm not sure what the standard bible would be nowadays for Signal Analysis books.

Reccomendation for M. Hamed, get pan (news client and email client) and
use eternal-september for free posting. I think that they do carry
alt.binaries.schematic.electronic (abse or a.b.s.e in common parlance).

Thanks. I'll check it out.

 

[josephkk]

I hardly ever use scilab, preferring octave, mainly because it makes a
better fist of running matlab stuff.
What I asked for wasn't graphics, but scilab code, ie. text.
OK.

No thanks, I climbed off the upgrade bandwagon aeons ago. Too many
dependency changes.

So not the point i wanted to make, that site is a script and cookie laden
hell. I quit instead of getting a download.

 

[josephkk]

ooOn Sat, 4 May 2013 13:42:51 -0700 (PDT), "M. Hamed"

I'm currently trying to go through "The scientist and engineer's guide to dsp". Most of these DSP books try to include very little mathematics and I feel this takes away from my understanding. I am thinking I should go for a Signals and Systems book for a better groundwork then move to DSP (except that it takes a long time this way ).

And you are correct on the low math aspect. Perhaps taking a senior or
graduate level course at a local uni or online would appropriate.

Back in my school days we had a Signal Analysis course and the book was Signal Analysis by Papoulis. It was a bit difficult to follow back then but it would probably be easier now. Can't find a good copy though.

Search engines are your friend.
Is this the book you mean:?
< http://www.amazon.com/Signal-Analysis-Athanasios-Papoulis/dp/0070484600

 

[Robert Macy]

Thanks for the tips on windowing. I'll take look into it.



You are absolutely right. I will look into some other alternative that shares files as files.

Here is the code:

///////////////////////
// Program Constants

start_freq = 550e3;     // starting frquency = 550 kHz
end_freq = 1700e3;      // end frequency   = 1700 kHz
freq_spacing = 20e3;    // channel spacing = 20 kHz
t_end = 1e-2;           // last time sample = 10 us
amplitude = 1;
sample_rate = 4e6;
random_phase = 1;
fftplot = 1;
nf = 20000             // number of frequency samples to plot

///////////////////////

// create time variable

time = 0:1/sample_rate:t_end;
szt = size(time, '*');

// create station frequencies spaced by 20khz
frq = start_freq:freq_spacing:end_freq;
// transpose
frq = frq';
szf = size(frq, '*');

//phase = (2*%pi/szf)2*%pi/szf):2*%pi;
phase = 2*%pi*rand(1,szf);
phase = phase'
phase_matrix = random_phase*phase*ones(1,szt);

// time frequency matrix, one row for each frequency and time is in columns
time_freq = amplitude*sin((2*%pi*frq*time)+phase_matrix);

// sum all the sinusoids
composite = sum(time_freq,;

// Now calculate and plot FFT
if fftplot then
    N = szt;

    y = fft(composite);

    f=sample_rate*(0N/2))/N;
    n=size(f,'*');

    clf()
    if nf=0 then
        nf = n
    end
    plot(f(1:nf),abs(y(1:nf)))
end

//////////////////////////////////////////////////

I don't use Scilab [I think that was its name] after a giant battle
with them over accuracy. They essentially told me tough.

After floundering for a few years, I discovered octave and have NOT
changed again. I do EXTENSIVE DSP work and octave not only predicted,
but taught me, a great deal about what was going on. During DSP
development after trying a process of manipulation in octave, it was a
piece of cake to write C/C++ real-time code to do the same thing. As
you know, battling code errors is enough to take on, not trying to
design the DSP processing at the same time.

Either you have to make certain you have a complete set of cycles for
each waveform in your FFT window, else you 'scatter' the energy of non
matching sinusoid into adjacent segments. Like 30% in one and 70% in
the next one, which makes your display look ghastly - but it is
correct by mathematical definitions. So...if you make the FFT window
of time match ALL waves THEN you'll get your expected results. If
impossible to match, you can use a window. Think of the window as
gently sliding into the data and then gently sliding back out, so the
'accidental' abrupt start and stop of your window over the waveforms
just doesn't happen. There are two very effective windows to use. go
to TI's ap note on their incredible ADS1282 24 bit data acquisition
system. They use a window that is a great compromise between
preserving energy and injecting noise into adjacent segments. In the
interim, I simply use an 'adjusted' hanning window which is fairly
easy and doesn't 'splatter' the energy around. Note, I said Hanning,
not Hamming. The shape of the hanning window is one cycle of a cos
wave, starts at zero goes to amplitude of 2, then returns to zero and
is EXACTLY the length of your FFT window. Octave provides that as a
supplied function so it's easy to use, but you can construct one
mathematically for your needs. Just keep in mind that whatever window
you use you MUST make the area under that window equal to one, else
you'll shift the amplitudes of your FFT, losing quantified answers,
only getting qualified answers.

Be FAR easier to talk in octave, rather than Scilab.