A digital demodulator for commercial FM broadcast

Material Information

A digital demodulator for commercial FM broadcast
Stevens, Jeffrey Donald
Publication Date:
Physical Description:
ix, 130 leaves : illustrations ; 29 cm

Thesis/Dissertation Information

Master's ( Master of Science)
Degree Grantor:
University of Colorado Denver
Degree Divisions:
Department of Electrical Engineering, CU Denver
Degree Disciplines:
Electrical Engineering
Committee Chair:
Ross, Douglas A.
Committee Members:
Anderson, Marvin F.
Thomas, Joe E.


Subjects / Keywords:
Digital communications -- Simulation methods ( lcsh )
Demodulation (Electronics) -- Simulation methods ( lcsh )
Radio -- Receivers and reception ( lcsh )
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )


Includes bibliographical references (leaf 87).
General Note:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Electrical Engineering.
Statement of Responsibility:
by Jeffrey Donald Stevens.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
22927287 ( OCLC )
LD1190.E54 1990m .S73 ( lcc )

Full Text
Jeffrey Donald Stevens
B.S., University of Wisconsin, 1985
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Master of Science
Department of Electrical Engineering

This thesis for the Master of Science degree by
Jeffrey Donald Stevens
has been approved for the
Department of
Electrical Engineering
Marvin F. Anderson
Date 2*7

Stevens, Jeffrey Donald
A Digital Demodulator for Commercial FM Broadcast
Thesis directed by Professor Douglas A. Ross
The world of commercial audio technology is
moving inexorably into the digital domain. Although
digital audio technology was in its infancy just a
decade ago, nowadays consumers can chose from dozens of
brands of digital Compact Disc players, Digital Audio
Tape recorders, Digital Signal Processors, and at least
one brand of digital preamplifier. To date, however,
there are no readily available digital FM radio
Most commercial FM radio receivers may be
conveniently divided into three circuit blocks, a tuning
block that has as its output a single radio channel at
an intermediate frequency (IF), a demodulator block that
recovers the audio information from the IF signal, and a
decoder block that separates the audio signal into its
left and right stereo channels. This paper emphasizes
the design and simulation of a digital implementation of
the demodulator block.
A brief history of digital audio technology is
presented along with the motivation for designing a
digital FM radio receiver. Some of the design tradeoffs
examined during the development of the digital demodula-
tor design are addressed. The digital demodulator

design and simulation programs are described and many
graphs and tables quantifying simulated digital demodul-
ator performance are analyzed in detail. Finally, con-
clusions regarding the feasibility of the digital demod-
ulator d,esign are presented.
The form and content of this abstract are approved. I
recommend its publicatic
as A. Ross

This thesis is dedicated to my wife, Pam, whose incredi-
ble patience and strength of character sustained me
throughout its generation, and to Steve Norris-Saucedo,
who may not realize that his friendship enabled me to
reach this point in my education.

I. INTRODUCTION ................................... 1
III. DIGITAL MODULATOR DESIGN ...................... 11
VI. CONCLUSIONS ................................. 85
REFERENCES ........................................ 87
B. FORTRAN PROGRAM LISTINGS .................... 91
C. PASCAL PROGRAM LISTINGS ..................... 122

1. Variables Used in Demodulator Block Diagram
Equations ...................................... 11
2. Cosine Integration Method Error Values ......... 50
3. Modem Simulation Initial Time Domain SNR
Values Double Precision ..................... 71
4. Modem Simulation Initial Time Domain SNR
Values Single Precision ...................... 75
5. Modem Simulation Initial Time Domain SNR
Values Random Phase, Single Precision ........ 77
6. Modem Simulation Initial Time Domain SNR
Values Single Precision, 0.01 Peak Noise .... 84
7. Modem Simulation Initial Time Domain SNR
Values Single Precision, 0.1 Peak Noise ..... 84

1. Block Diagram of FM Radio Receiver Showing
Digital Demodulator Stage ..................... 8
2. Block Diagram of FM Radio Receiver
Demodulator Stage ............................. 12
3. Block Diagram of FM Radio Receiver
Demodulator Stage Digital Implementation .... 18
4. FM Spectrum of a Single Baseband Frequency .... 23
5. FM Spectrum of a Single Baseband Frequency
Showing Carrier Suppression .................... 25
6. FM Spectrum of Two Baseband Frequencies ........ 26
7. FM Spectrum of Five Baseband Frequencies
with Arbitrary Phase Angles ................... 27
8. FIR Lowpass Filter Response Blackman
Window Function ................................ 31
9. FIR Lowpass Filter Response Chebyshev
Window Function ................................ 35
10. FIR Filter Phase Delay Verification ........... 42
11. Modem Simulation Phase Delay Verification ..... 46
12. Modem Simulation Output Spectrum Example ...... 48
13. FM Spectrum for Fs = 400kHz, Fc = 100kHz,
Fd = 75kHz, and Fm = 15kHz .................... 53
14. FM Spectrum for Fs = 800kHz, Fc = 200kHz,
Fd = 75kHz, and Fm = 15kHz .................... 54
15. FM Spectrum for Fs = 802816Hz, Fc = 200704Hz,
Fd = 75kHz, and Fm = 12544Hz .................. 56

16. FM Spectrum for Fs = 999424Hz, Fc = 249856Hz,
Fd = 75kHz, and Fm = 15616HZ .................. 58
17. FM Spectrum for Fs = 999424Hz, Fc = 249856Hz,
Fd = 50kHz, and Fm = 15616Hz .................. 60
18. Revised I and Q Channel FIR Lowpass Filter
Frequency Response ............................. 61
19. Arctangent Block Output Time Series ............ 63
20. Modem Simulation Time Domain Output with
Differentiator Delay Value Too Large ........... 64
21. Modem Simulation Time Domain Output with
Differentiator Delay Value Acceptable .......... 66
22. Modem Simulation Spectrum without Pre-
decimation Filter.................................... 67
23. Frequency Response of Pre-decimation FIR
Filter Chebyshev Window Function ............. 69
24. Modem Simulation Spectrum with Pre-
decimation Filter.................................... 70
25. Modem Simulation Spectrum for Seven
Baseband Harmonics ............................. 72
26. Modem Simulation Time Domain Output for
Seven Baseband Harmonics ....................... 74
27. Modem Simulation Time Domain Output for Seven
Harmonics, Random Phase, Single Precision ...... 78
28. Modem Simulation Spectrum for Seven Harmonics,
Random Phase, Single Precision ................. 79
29. FM Time Series Spectrum for Seven Harmonics,
Random Phase, Random Noise ..................... 81
30. Modem Simulation Spectrum for Seven Harmonics,
Random Phase, Single Precision, Random Noise .. 83

Why design an FM radio receiver utilizing digi-
tal technology when analog radios offering good perfor-
mance at low cost are readily available? The motivation
can be partially derived from the following cursory
examination of the recent history of commercial audio.
The conversion of American home audio systems
from analog to digital technology began in the early
1980s with the introduction of a Compact Disc (CD)
player manufactured by Sony. Initially, public accep-
tance of the new medium was slowed by the limited selec-
tion of prerecorded music. In addition, many music
critics were quick to denounce the "digital sound" of
the first generation of CD players as harsh, grating,
and generally unlistenable. Prices of $1000 and more
for these early players prevented a large segment of
audio consumers from becoming familiar with the new
technology. Despite these and other obstacles, CD
players gained a solid foothold in this country.
As the recording industry and component manufac-
turers climbed the digital audio learning curve, fewer

critics complained about the sound quality of CDs.
Improved yields from disc producers allowed lower prices
and an expanded catalog of artists and titles. CDs
began to crowd analog LP records from the shelves of
music shops across the country. CD player prices
dropped sharply from $1000 and up to less than $500 for
many of the second generation models. The digital con-
version initiated by the CD player was well underway.
The second commercial audio product to embrace
digital technology was the Digital Audio Tape recorder
(DAT). Political pressure has delayed the introduction
of DAT recorders in this country, but they are widely
available in Japan, and can be bought in America from
"gray market" importers. DAT recorders combine the same
superb sound quality that can be heard from CD players
with the capability of home digital recording. A home
digital recording of a CD or DAT format selection should
result in an essentially perfect copy of the original
material. Once all of the copyright issues are worked
out, DAT recorders will probably become as common as CD
players are now.
With recorded music mediums taken care of, audio
equipment manufacturers turned their attention to the
powerful signal processing capabilities available from
digital logic technology. The next wave of digital
audio consumer products were the Digital Signal Proces-

sors (DSP), which allowed home listeners to simulate the
acoustic properties of live music venues ranging from
massive cathedrals to intimate clubs. As was the case
with the CD and DAT machines before them, prices for the
DSP products are dropping rapidly even as the units are
becoming more sophisticated.
The most recent entrant into the digital audio
field is the digital preamplifier. These devices use
digital circuits for volume, balance, and tone control
as well as for switching among digital sound sources. A
system using a digital preamplifier in concert with the
equipment described above would allow the music repro-
duction process to be done in the digital domain up to
the point of power amplification. In other words, the
sole stereo pair of digital-to-analog converters for
such a system could be placed in the power amplifier
The only widely used high fidelity music source
missing from the above audio system is an FM radio
receiver. Therefore, it is virtually inevitable that
digital FM radio receivers will become commonplace.
Of course, inevitability in itself is not neces-
sarily the most intellectually satisfying reason for
designing a new product. Digital design techniques
offer several inherent advantages over their analog
counterparts, including better noise immunity, freedom

from component performance drift due to aging and
temperature effects, the linear phase response from
finite impulse response (FIR) digital filters, and the
availability of sophisticated, yet inexpensive signal
processing integrated circuits.
In addition, the digital FM radio receiver
offers the prospect of adaptive digital filtering for
noise suppression, digital signal processing for multi-
path cancellation, and fast switching among several
diversity antennae in mobile applications to name just a
few potential design opportunities. It seems clear that
the digital FM radio receiver is not only inevitable,
but highly desirable as well. This paper focuses on a
crucial first step towards the realization of a digital
FM radio receiver: the design and simulation of a high
fidelity digital FM demodulator.
The discussion begins with an overview of some
of the design tradeoffs considered in the early stages
of the design process. The chosen design is examined at
the block diagram level of detail. The development of
baseband and IF simulation programs for both time and
frequency domain analysis of radio performance is
presented. The effects of calculation precision and
receiver noise on radio performance are investigated.
Finally, conclusions are reached on the feasibility of
implementing the digital demodulator in hardware.

The first question faced by the digital radio
designer can be stated in a deceptively simple form:
"Where should I place the analog-to-digital converter
(A/D)?" Theoretically, at least, the A/D can be placed
anywhere in the radio receiver except ahead of the
receiving antenna. Indeed, the options are so numerous
and the factors involved so extensive that an exhaustive
study of them would form a thesis topic in itself. Only
a few of the possibilities will be discussed here.
The most desirable location for the A/D in the
eyes of a digital designer is as close to the antenna as
possible. This location results in the replacement of
the greatest number of analog circuit elements with
their digital counterparts. In this case, referred to
here as RF sampling, the tuning and frequency downcon-
version processes as well as signal demodulation would
be performed in the digital domain. While these pro-
cesses would not be prohibitively difficult to attain
with commercially available integrated circuits (ICs),
the exceptionally high conversion rate required of the

A/D is simply not widely available in today's main
stream technologies.
The commercial FM radio band extends from 88
megaHertz (MHz) to 108MHz (roughly). To satisfy the
Nyquist sampling criterion, the A/D would have to gener-
ate samples at a rate of at least 216MHz. Real world
hardware requires a safety margin over the Nyquist rate,
which means that the actual minimum sample rate is
probably somewhere in the vicinity of 250MHz. While
there are a few A/D converters that operate at and above
this rate, they are very expensive and limited to 8 bits
of resolution at best. The examination of whether or
not this is sufficient resolution for a high fidelity FM
receiver is somewhat beyond the scope of this paper. In
any case, the high cost of the A/Ds was deemed suffi-
cient to rule out this approach.
A second approach to RF sampling was also brief-
ly considered. If the 20MHz wide FM band were bandpass
filtered prior to the A/D conversion, the bandpass
sampling technique could be utilized. This technique
requires a conversion rate that is at least two times
the input bandwidth. For the FM band, a conversion rate
of 56 MHz would satisfy the bandpass sampling criterion
and also place the copies of the FM spectrum exactly
halfway between multiples of half the sample rate
(Fs/2). Before any significant effort was expended on

behalf of this approach, it was also ruled out due to
excessive A/D cost.
Baseband sampling lies at the opposite end of
the conversion frequency spectrum (as well as the
sophistication spectrum) from the methods described
above. In this somewhat trivial case, the audio output
of a conventional FM radio receiver is sampled by an A/D
converter. While this design could incorporate a digi-
tal stereo decoder, it breaks little new ground and was
therefore ruled out due to lack of interest.
With RF and baseband sampling both out of the
running, the only place left to locate the A/D was some-
where in the IF stage of the radio. Since it was
desired to maximize the use of digital techniques in the
demodulator, the A/D was located at the head of the IF
circuit block. A block diagram of the proposed overall
receiver concept is presented in Figure 1.
The workings of the analog input amplifier,
tuner/downconverter, low pass filter and automatic gain
control (AGC) blocks will not be described here. What
follows is a brief explanation for their inclusion. The
input amplifier is needed to boost the weak antenna
voltage signal (typically a few microvolts to a few
millivolts) to levels compatible with the requirements
of the tuner block. The tuner/downconverter isolates a
single station from the FM band and translates it down

Figure 1
Block Diagram of FM Radio Receiver Showing
Digital Demodulation Stage

to IF. The filter is needed to bandlimit the highest
signal passed to the A/D in deference to the Nyquist
sampling theorem. Finally, the AGC is incorporated to
make full use of the A/D's available resolution. This
minimizes quantization noise and also prevents unde-
sirable volume variations when changing radio stations.
Conventional analog FM receivers utilize a stan-
dard IF of 10.7 MHz, a frequency for which all of the
analog modulation circuit elements can be optimized.
Digital components, unlike their analog counterparts, do
not require optimization to any particular frequency
band. Therefore, the IF for the radio shown in Figure 1
can be decreased from 10.7 MHz to the minimum IF needed
to support the bandwidth of a single channel.
Assuming a channel bandwidth of 200kHz plus a
guard band 50kHz, a simulated IF in the neighborhood of
150kHz was deemed appropriate. Compatibility with the
established CD player sample rate of 44.1kHz was con-
sidered to be an essential criterion in the selection of
the digital radio sample rate, which in turn influences
the selection of the IF. To make the most effective use
of the guard bands, it was desired to have the IF signal
centered between 0Hz and (Fs/2)Hz.
An IF of 176.4kHz with a sample rate of 705.6kHz
meets all of the above requirements. The sample rate is
sixteen times the CD sample rate of 44.1 kHz, the IF is

centered between 0Hz and 352.8kHz (Fs/2), and the guard
bands are 76.4kHz wide.
Although the block diagram shows an FM stereo
decoder implemented in the digital domain, time con-
straints led to the decision not to include its design
in this thesis.

A block diagram of the demodulator stage is pre
sented in Figure 2. Mathematical equations are used to
describe the operation of each of the demodulator
blocks. Table 1 is a list of the variables used in the
X(t) = demodulator input
A = maximum amplitude of demodulator input
w = IF in radians/second
0(t) = IF phase angle in radians
0'(t) = time derivative of the IF phase angle
Y(t) = demodulator output
Table 1
Variables Used in Demodulator Block
Diagram Equations
The following equations are presented to demon-
strate the function of each of the demodulator blocks in
mathematical terms. The input to the demodulator
circuit is given by
X(t) = ACos(wt + 0(t)). (a)
The outputs of the I and Q channel multipliers are
I(t) = A[cos(wt + 0(t))sin(wt)]
Q(t) = A[cos(wt + 0(t))cos(wt)].

Input '
Figure 2
Block Diagram of FM Radio Receiver
Demodulator Stage

Using trigonometric identities, these two equations can
be rewritten to call out the sum and difference frequen-
cies explicitly:
I(t) = 3*A[sin(2wt + 0(t)) sin(0(t))] (c)
Q(t) = ^A[cos(2wt + 0(t)) + cos(0(t))].
The I and Q channel low pass filters (LPFs) remove the
sum frequency components located at two times IF which
results in
ILPF(t) = -%Asin(0(t)) (d)
QLPF(t) = liAcos(0(t)) .
The output of the ratio generator block is
R(t) = -tan(0(t)). (e)
The output of the arctangent generator block is
Arc(t) = -0(t). (f)
The output of the differentiater block is
Y(t) = 0'(t), (g)
which is seen to be the recovered original modulating
signal. (Appendix A contains a brief derivation of the
equation governing FM transmission which demonstrates
that the baseband audio signal is contained in the
derivative of the carrier phase angle.)
Recalling that the demodulator operates in the
digital domain, it is clear that all of the circuit
blocks must have digital implementations. The following
paragraphs describe the particular methods selected to
create each of the demodulator circuit blocks.

The I and Q channel multipliers each simply
multiply their two binary inputs together to produce
single binary outputs according to the equation
Y = A B,
where A and B are the inputs and Y is the output.
Unlike their analog counterparts, the digital multi-
pliers are linear devices, which means that their out-
puts contain no terms of the type
Y = nA mB,
where n and m are integers. Thus, the multipliers
generate only the Sum and difference frequencies shown
in equation (c) on the previous page with no added
spurious components. This is one of the advantages
inherent in digital design techniques.
The sum frequencies from the I and Q channel
multipliers shown in equation (c) carry redundant infor-
mation and are, therefore, filtered out. Linear phase
finite impulse response (FIR) low pass digital filters
are employed for this purpose. The transfer functions
of the FIR filters are chosen to pass the difference
frequency shown in equation (c) and to attenuate the sum
frequency. Identical filters are used in the I and Q
channels of the demodulator. The ability to generate
steep transition band filters with linear phase charac-
teristics is a second advantage enjoyed by the digital
demodulator over analog designs.

The ratio generator simply calculates the ratio
of the I and Q channel output signals to produce the
tangent of the IF phase angle. As with the preceding
circuit elements, the ratio generator is a linear device
which introduces no spurious frequency components. Note
that the generator will require a means to avoid divi-
sion by zero. Division by a variable identically equal
to zero is highly unlikely in a computer simulation, so
the potential was ignored in the discussion of the simu-
lation program presented in the next chapter. However,
the threat is greater in a hardware implementation due
to the limited quantization levels available. One
solution would be to always place the larger value in
the denominator of the ratio generator. The output
could then be inverted as necessary.
The arctangent generator calculates the arctan-
gent of its input. Since the arctangent operation is
not a single valued function, that is, it could generate
an infinite number of correct outputs for a given input,
care must be exercised in its implementation. Specifi-
cally, to avoid output angle ambiguities, the sample
rate must be chosen to ensure that the arctangent gener-
ator output could not increment by a value larger than tt
from one sample period to the next. In addition, care
must be taken to avoid erroneous results when the
arctangent moves between quadrants II and III in succes-

sive sample periods. These two concerns will be
explored more thoroughly in the next chapter.
Because the demodulator operates in the digital
domain, the differentiation block can only be approxi-
mated by a finite difference method. In this method,
the current arctangent block output value is subtracted
from a previous arctangent output value that has been
delayed by one or more sample periods. Division by a
small increment of time is inherent in sampling process,
so the incremental change,in arctangent (delta vided by the incremental change in time (delta t) ap-
proximates the time derivative of the arctangent. Note
that the expected order of subtraction is reversed in
the above description to cancel the negative output of
the ratio generator. (Although the need to avoid a 180
degree input to output phase shifts in audio equipment
has not yet been resolved to everyone's satisfaction,
effort was made to preserve absolute phase in the dig-
ital demodulator.)
The output filter is needed to block the ultra-
sonic frequency components that are a direct result of
the sampling process. While these frequencies are
inaudible, they can interfere with tape recorder bias
frequencies and overload downstream analog amplifiers.
As was the case for the I and Q channels, an FIR digital
filter was chosen for its linear phase characteristics.

The coefficients were chosen to produce a flat response
out to 15kHz and about 40dB of attenuation at the FM
pilot tone frequency of 19kHz.
Not shown in Figure 2 is a decimation circuit
block used to divide the digital demodulator output
sample frequency down to 44.1kHz, the rate used in CD
players. The decimation block isn't shown because it is
not required for demodulator operation and because the
block diagram was meant to lend itself to analog as well
as digital implementation and the decimation block has
no analog counterpart.
A summary of the digital blocks used to imple-
ment the FM demodulator is shown in Figure 3.

Figure 3
Block Diagram of FM Radio Receiver Demodulator
Stage Digital Implementation

Several programs were written during the course
of this research to simulate the time and frequency
domain characteristics of signals at baseband and at IF,
as well as to generate coefficients for and characterize
the response of digital FIR filters. All of the simul-
ation programs and the FIR filter programs were written
in Fortran because of its wide selection of intrinsic
functions and its complex number capabilities. All of
the programs were written to be interactive with prompts
asking the user to enter values for frequencies,
amplitudes, filter lengths, window selection, and the
like. A listing of Fortran source code for a represen-
tative of each type of program used to research this
thesis can be found in Appendix B.
One problem with Fortran is the lack of flexible
graphics functions. Since graphical output was of fun-
damental importance in the debugging of simulation
programs and the selection of various demodulator para-
meters, it was decided to write the Fortran program
output data to ASCII files that could be read into

plotting routines written in another language. The
specific language selected was Turbo Pascal. Coupled
with its associated Graphics Toolbox collection of
routines, Turbo Pascal provided the capability to create
the accurate, thoroughly annotated graphs presented in
this paper. The source code for a sampling of the
graphics programs is listed in Appendix C.
The objective of this phase of the thesis was to
write a program that simulated the operation of the
digital demodulator as shown in Figure 3. But before
the digital demodulator operation could be simulated,
data representing a frequency modulated baseband signal
had to be generated. This data, which would become the
input to the demodulator simulation program, could be
created with a program based on a digital interpretation
of the equation defining frequency modulation. (A deri-
vation of such an equation is presented in Appendix A.)
The first program to be written, called
FM1C0S.FOR, simulated the frequency modulation of a
single cosine with variable frequency and amplitude and
a 0 phase angle onto a variable frequency carrier. An
embedded FFT routine allowed the FM spectrum of any
combination of modulation frequency and amplitude, car-
rier frequency, and sample frequency to be graphed.
FM1COS.FOR is made up of four basic blocks. The
first of these blocks declares and initializes the vari-

ables used in the program, either by direct assignment
or via interactive keyboard data entry. Two main cal-
culation loops follow the variable initialization block.
The first of these loops calculates the values of an N-
point time series representing the amplitude of a fre-
quency modulated baseband cosine uniformly sampled at a
rate given by Fs. The value of N can be varied by the
user, but the FFT routine requires N to be equal to an
integer power of two.
The last continuous time equation shown in the
derivation in Appendix A involves the integral of the
baseband modulating signal. When the modulating signal
is a single cosine, the value of the integral can be
calculated in one of two ways. The integral can be
evaluated to give a sine function which is then used in
the time series calculation. Alternatively, the
integral can be approximated with a discrete time sum-
mation of cosine terms.
The FM1C0S.F0R program is set up to calculate
the amplitude time series values using both methods.
Because it is mathematically exact, the time series
value based on the sine function method is considered to
be the correct value at each time step. An error term,
equal to the average squared difference between the two
methods, is calculated at each time step. This error
term, expressed on a decibel scale, is written to the

screen. Error values corresponding to various
combinations of baseband cosine and maximum deviation
frequency are reported in the next chapter.
The second loop in FM1C0S.F0R calls an FFT
subroutine (described below) and then uses the subrou-
tine output to calculate the spectral power at each
frequency step. The N/2+1 spectral power values are
written to an ASCII file that is used to generate plots
of frequency versus power. Along with the spectral
power values, FM1C0S.F0R also writes the sample, base-
band and carrier frequencies and the length of the time
series to the ASCII file for graph annotation.
The final element in FM1C0S.F0R is the FFT
subroutine, which was taken verbatim from reference [1].
The subroutine, called SPFFTC, uses a symmetric Cooley-
Tukey algorithm based on time decomposition with input
bit reversal.
Once written, FM1C0S.F0R needed to be verified.
By careful manipulation of baseband frequency and ampli-
tude, carrier frequency, and sample frequency values, FM
spectra published in the communication theory text by
Carlson could be duplicated [2]. An example showing an
FM1C0S.FOR output that closely corresponded to one of
the graphs presented in reference [2] is shown in Figure
4. The amplitude scale in Figure 4 is linear to match
the scale used in the reference.

flHplitudt 05. Frequency
Figure 4
FM Spectrum of a Single Baseband Frequency

A second method of verification was to check for
complete carrier suppression for beta values of 2.4 and
5.5. Figure 5 is a graph demonstrating complete carrier
suppression for a beta value of 2.4. Together, these
two verification methods confirmed the correct operation
of the modulation simulation program.
With the FM1C0S.FOR simulation program working
correctly, it was a simple matter to add a second cosine
with variable frequency and amplitude and a 0 phase
angle to form the program called FM2C0S.F0R. As with
FM1C0S.FOR, this new program was verified by comparing
its output with graphs found in reference [2].. An
example of FM2C0S.F0R output plotted against a linear
amplitude scale is shown in Figure 6.
The final refinement added to the modulator
simulation program before moving on to demodulator
simulation was a provision for frequency modulating any
number of baseband cosines, each with pseudorandom phase
relative to the others. In this program, FMNCOS.FOR,
each phase term was generated by using the pseudorandom
number routine, SPRAND, listed in reference [1]. This
subroutine outputs pseudorandom values between zero and
one. The phase values were derived by subtracting 0.5
from the SPRAND output value and multiplying the result
by 27T. Figure 7 is an FFT generated by FMNCOS.FOR with
an input of 5 baseband cosines.

flHHitude us. Frtautnty
Figure 5
FM Spectrum of a Single Baseband Frequency
Showing Carrier Suppression

Hup limit vs. Frtdutficy
Figure 6
FM Spectrum of Two Baseband Frequencies

Figure 7
FM Spectrum of Five Baseband Frequencies
with Arbitrary Phase Angles

With flexible, functioning modulator programs
ready to generate input data for the demodulation simu-
lator, it was time to begin writing digital demodulator
simulation code. In order to achieve the highest pos-
sible fidelity, it was decided to write the demodulator
code using double precision variables and functions.
Once the demodulator concept was proven valid, the code
could be backed down to single precision. To support
the double precision demodulator program, FM1C0S.F0R,
FM2C0S.F0R, and FMNCOS.FOR were converted to double
precision as well. The three single precision precur-
sors were retained, and the six programs were renamed
FMxCOSSP.FOR or FMxCOSDP.FOR (x = 1, 2, or N) to avoid
In addition to the conversion to double preci-
sion, one other enhancement was added to the modulator
simulation programs. Part of the digital demodulator
research agenda was to investigate its performance given
a noisy input signal. With that in mind, the modulator
programs were modified to provide noise in the form of a
variable amplitude pseudorandom signal (PRS) added to
the IF time series.
The pseudorandom sequence is supplied by the
SPRAND function found in reference [1]. A zero mean
sequence is created by subtracting 0.5 from every SPRAND
output (X). The maximum desired amplitude (A) attain-

able by the PRS is entered by the user in response to a
program prompt. The added noise signal is given by
Noise = A(X 0.5). The noise term can be eliminated by
entering zero when prompted for the maximum noise level,
but this will result in a nonfatal runtime error.
Rather than write the modulator simulation
program output to an ASCII file that would then be read
into the demodulator program, it was decided to incor-
porate both modulator and demodulator simulation code
into a single program. As was the case for the series
of modulator programs, the modulator/demodulator (Modem)
programs were written to accommodate first one, then
two, and ultimately an arbitrary number of baseband
cosines. In the discussion of the demodulator simu-
lation code that follows, the code segments will be
covered in the same order as the circuit blocks (as
shown in Figure 3) they represent.
Discounting the generation of discrete time
samples via A/D conversion, which is a process inherent
in computer simulations, the first simulated circuit
elements in the digital demodulator are the I and Q
channel multipliers. There is little to say about them,
however, since they can each be simulated by a single
line of code that multiplies the modulator block output
value at each time step by the value of either the sine
or cosine of the carrier frequency, respectively, at

that time step.
The I and Q channel FIR low pass digital filters
required more effort. Work on the Modem programs, per
se, was postponed pending the creation of an FIR design
program. Reference [1] provided the kernel, called
SPFIRD, of a sync-function-based single pass FIR filter
design routine with the option of choosing a rectan-
gular, triangular, Hamming, Hanning, or Blackman window.
The SPFIRD subroutine provides for the design of low-
pass, highpass, bandpass, or bandstop filters.
A main program was written that prompts the user
for the sample frequency, cutoff frequency or frequen-
cies, filter type, order and length, and window function
selection. After the filter coefficients are calculated
by SPFIRD, another subroutine from reference [1], called
SPGAIN, is used to find the complex gain of the filter.
The filter's frequency response is calculated from the
complex output of SPGAIN. A plotting program was
written to display filter frequency response along with
some of the user-entered parameters. The calculated
filter coefficients and filter length are written to one
of three ASCII files for use by the Modem programs.
Figure 8 is a graph of the transfer function of
one of the first FIR filters designed for the I and Q
channels of the demodulator. Note that filter atten-
uation is plotted on a logarithmic scale because it

ahp I itudt 'is Frequency ~
Figure 8
FIR Lowpass Filter Frequency Response
Blackman Window Function

provides better resolution of the filter's stopband
attenuation than does a linear scale. The 75dB of
attenuation afforded by this filter, which uses a
Blackman window, is the greatest amount of attenuation
that can be achieved with a program based on the kernel
supplied by reference [1].
It was decided to perform all of the Modem
program filtering operations in the time domain using
loops to directly implement the convolution summation
equation defining an FIR filter. Because of the need to
avoid the transient response phase of the filter pro-
cess, the filter input time series had to be lengthened
by a number of samples equal to the order of the filter.
To elaborate, if the filter comprised M coefficients and
the desired filter output time series was to be N
samples long, the length of the filter input time series
is given by (N+M-l), where (M-l) is the order of the
The first element in the filter output time
series is equal to the convolution summation of the
first M elements of the filter input time series with
the M filter coefficients. Subsequent output time
series elements are calculated in an analogous manner.
This method of avoiding filter transient response by
extending the filter input time series is used through-
out the Modem programs. An example of the method for an

FIR filter of order two is shown below. The variable X
denotes filter input time series elements, B denotes
filter coefficients, and Y denotes filter output time
series elements. Only the first two output series
element calculations are shown. Note that the selection
of time series indices is completely arbitrary. This
method accurately simulates steady state demodulator
operation, which would be achieved within much less than
one second following circuit turn-on.
X(-2) X(-l) X(0) X(l) X(2) ...
B(0) B(l) B(2)
Y(0) = X(0)B(2)+X(-1)B(1)+X(-2)B(0)
X(-2) X(-l) X(0) X(l) X(2) ...
B(0) B(l) B(2)
Y(1) = X(1)B(2)+X(0)B(1)+X(-1) B(0)
During the process of improving the performance
of the first Modem program, it was determined that 75dB
of attenuation was probably not sufficient to permit
optimum performance of the demodulator as simulated with
double precision variables and functions. Research into
window function characteristics led to an FIR filter
design program called WINDOW written by Rabiner, McCone-
gal, and Paul that provided the window functions listed
above as well as Kaiser and Chebyshev window functions.
The Fortran source code for WINDOW and a brief explana-

tion of its operation can be found in reference [3].
The modified WINDOW source code as used for this thesis
is listed in Appendix A.
Use of the Chebyshev window allows for two of
the following three parameters to be independently
specified: amplitude ripple, filter length, and nor-
malized transition band width. The non-specified third
parameter is calculated by the program. Selection of
the Chebyshev window resulted in the achievement of
stopband attenuations on the order of 120dB to 130dB.
The effects of the improved attenuation will be dis-
cussed in the next chapter. Figure 9 is a graph of the
frequency response of an FIR filter using the Chebyshev
window function. Note that the filter length and cutoff
frequency are the same as those for the filter whose
response is shown in Figure 8.
Up to this point in the discussion, the I and Q
channel signals have been multiplied by the sine and
cosine of the carrier frequency, respectively, and the
resulting sum frequencies have been filtered out. The
next step in the demodulator is to form the ratio of the
remaining I and Q channel difference frequencies. This
ratio is the tangent of the time varying carrier phase
term. Next, the arctangent of the ratio block output is
calculated. This calculation isolates the time varying
carrier phase term which contains the audio baseband

HMiituat vs. Frjoutficy
Figure 9
FIR Lowpass Filter Frequency Response
Chebyshev Window Function

information as shown in the derivation in Appendix A.
The simulation of the two demodulator steps
described above requires only a single line of code.
The Fortran intrinsic function ATAN2(argl,arg2) calcu-
lates the arctangent of the ratio of its two arguments,
(argl/arg2), automatically placing the resultant angle
in the correct quadrant based on the signs of the argu-
ments. As noted previously, since the arctangent opera-
tion is periodic, it is not a function; that is, a given
input does not have a unique output. Rather, there are
an infinite number of periodic output values associated
with every input value.
If managed carefully, the fact that the arctan-
gent is not a function doesn't cause a problem for the
digital demodulator. This is because the demodulator
relies on changes in phase of the IF carrier rather than
its absolute phase to recover the baseband audio infor-
mation. If the changes in phase are kept small and if
the quadrant II/quadrant III phase transitions are
accounted for, the digital demodulator will accurately
recreate the audio baseband signal.
An output signal corresponding to these changes
in IF signal phase is generated in the finite difference
segment of the digital demodulator simulation program.
The difference is calculated by subtracting the value of
the current arctangent block output sample from the

value of a previous arctangent block output sample that
has been stored. The storage time, or delay, is entered
by the user in response to a program prompt.
As was the case for the FIR filters, the opera-
tion of the finite difference block required lengthening
the input time series. Using the terminology from the
preceding paragraph, suppose the current arctangent
block output sample is the last sample in the arctangent
block output time series. Then the difference calcu-
lated by subtracting this last sample from a previous
arctangent block output sample becomes the last sample
in the finite difference block output time series. That
is, the last sample out of the finite difference block
has the same index as the last sample into the block.
In like manner, the first sample out of the
finite difference block has the same index as the first
sample into the block. The problem is this: when the
current arctangent block output sample is the first
sample in the arctangent block output time series, the
finite difference block has no previous value stored to
subtract it from. This previous value must be supplied
by lengthening the Modem baseband time series by a
number of samples equal to the largest expected finite
difference block delay value.
This delay, along with the sample rate, are the
two factors that must be carefully managed if the base-

band signal is to be accurately demodulated. Together
they determine whether or not the IF phase angle dif-
ference output by the finite difference block can exceed
2n. The selection of sample frequency to avoid aliasing
was discussed in a previous chapter. The selection of
an appropriate delay value was an empirical process that
will be discussed in the next chapter. At this point,
suffice it to say that the digital demodulator provided
a high fidelity output given suitable choices for sample
frequency and delay values.
There is an additional factor affecting the
finite difference block performance that must be ad-
dressed. If the difference between an angle in quadrant
II and another in quadrant III is to be calculated, the
result will be in error by a magnitude of 2ir unless a
correction term is added. This occurs because the
Fortran AT AN 2 intrinsic function assign values from ir/2
to 7T to angles in quadrant II and values from -ir/2 to -ir
to angles in quadrant III. If the current angle in the
differentiator block has a value of 0.9w and the delayed
angle has a value of -0.97T, the calculated output value
would be -0.9ir 0.97T, or -1.87T. The correct value is
seen to be 0.2ir. (Recall that the finite difference
block is "backwards" to cancel the negative output of
the tangent block.)

The digital demodulator simulation program
corrects for this calculation error by testing the value
of the finite difference block output. If the output
value is less than -it, the program assumes that the
phase is changing in a clockwise direction. In this
case, an example of which appears in the previous para-
graph, 2ir is added to the output. Similarly, if the
output value is greater than ir, the program assumes that
the phase is changing in a counterclockwise direction.
To correct this error, 2n is subtracted from the angle.
Finally, note that this quadrant factor imposes
a more stringent requirement on the delay variable than
does the fact that the arctangent is not a function.
The value of the delay must be chosen such that the
magnitude of the expected difference calculated by the
finite difference block is less than n. If the expected
difference magnitude exceeds tt, the correction calcu-
lations described above will corrupt it.
The output of the finite difference block is the
baseband audio signal represented as a time series of
samples. While in theory the digital demodulation
process is complete, the sample rate at this point,
initially chosen to be 705.6kHz, far exceeds what is re-
quired to satisfy the Nyquist sampling theorem. Indeed,
it is sixteen times the standard CD sample rate of
44.1kHz. The next two blocks in the digital demodulator

simulation program are included to decimate the sample
rate down to 44.1kHz.
Because decimation decreases the sample fre-
quency, spurious noise components that were above the
highest frequency in the audio band (and therefore
inaudible) but below Fs/2 may alias into the audio band.
To avoid this potential problem, a linear phase FIR
digital filter was placed ahead of the decimation block.
This filter was implemented in the time domain in the
same manner as the I and Q channel filters described
above. It was designed to have a frequency response
attenuated less than ldB at 15kHz and about 40dB at the
FM pilot tone frequency of 19kHz.
The decimation block is the final block in the
digital demodulator circuit. It is simulated by simply
writing every Nth sample to the demodulator output
array, where N is the decimation factor. The variable N
is entered by the user in response to a program prompt.
Because of the desire to perform an FFT on the demodu-
lator output, N is required to be a power of two.
Values of N less than sixteen are supported to allow the
display of any spurious frequency components out to the
undecimated Fs/2.
In order to assess the fidelity of the modula-
tor/demodulator (Modem) simulation it was necessary to
devise some means of comparing the demodulator baseband

output to the modulator input cosines. It was decided
to perform the comparison in the time domain by superim-
posing the input and output waveforms. This requires
both the artificial alignment of the two series in time
and the normalization of their amplitudes.
The first requirement can be satisfied by the
determination of the total phase delay caused by the
Modem processes. The only two processes that cause
phase delay in the Modem simulation are the FIR filters
and the finite difference block.
The phase delay resulting from passing data
through a linear phase FIR filter is known to be equal
to half the filter order. A program was written to
verify that the filters used in the Modem simulation
imposed the predicted phase delay. The program,
DELAYTST.FOR, implements an FIR filter in the time
domain, with input and output time series indexed as
described above for the I and Q channel filters. The
two time series are written to an ASCII file that is
read by a Pascal plotting routine written specifically
to display them. Figure 10 shows the ten sample delay
caused by an FIR filter of order twenty. Note that the
graph also reveals the expected 3dB of attenuation at
the cutoff frequency of the filter. This and other
examples verified the phase delay performance of the FIR
filters used in the Modem simulation.

Figure 10
FIR Filter Phase Delay Verification

The calculation of phase delay caused by the
finite difference process is analogous to the FIR filter
case. The phase delay is equal to half the order of the
finite difference delay, where the order is equal to one
less than the difference in indices of the two samples
being operated on.
The total Modem simulation phase delay is calcu-
lated by adding the phase delays of its components.
Because the I and Q channel FIR filters are in parallel
and are identical, their combined phase delay is equal
to their individual phase delays. Let the order of the
I and Q channel FIR filters be denoted 0,. Similarly,
let the orders of the pre-decimation FIR filter and the
finite difference block be denoted 02 and 03, respec-
tively. Then the total phase delay, d, of the Modem
simulation program is given by
= (01+02+03)/2.
This calculated total phase delay was used to offset the
index of the Modem program output time series in the
time domain plotting program.
The second step required to accurately superim-
pose the input and output time series is correct ampli-
tude normalization. Input normalization for the Modem
programs lacking pseudorandom baseband cosine phase
angles is achieved by simply dividing the input cosine
time series amplitude calculated at each point by the

number of cosines in the baseband signal. Output time
series normalization required more thought.
It was decided that there was no analytical way
to predict what the maximum output amplitude would be
for all possible demodulator operating parameter values.
Rather, all members of the output time series would be
divided by a normalization factor determined empirically
for each execution of the simulation program.
A least squares method was selected to generate
the required normalization factor. The factor is calcu-
lated by finding the ratio of two summations. The
numerator is a summation of the products of the input
time series value and the output time series value whose
index has been offset according to the discussion above.
The denominator is a summation of squared values of the
output time series. Both summations are calculated over
the entire range of output time series indices. The
output time series is normalized by multiplying each
member of the series by the normalization factor.
This normalization method offers two advantages
over a more common, simpler method, which is to search
the output time series for its maximum value and then
divide each member of the series by that value. The
first is its insensitivity to an isolated large time
series value which would distort the results of the
simpler method. The second is its ability to normalize

a low frequency output cosine. The simpler method
requires the presence of time series values representing
at least one half of a cosine period of lowest output
frequency. This is to ensure that a true maximum value
is available in the series. Since the selected method
compares the output to the input rather than searching
for a maximum, it can correctly normalize an output
frequency represented by a set of time series values
covering less than one half of the associated period.
With the above two factors in mind, time domain
plots were generated that superimposed the offset and
normalized output time series onto the normalized input
time series. If the calculated phase delay and nor-
malization are correct and the Modem simulation program
is operating correctly, the input series and offset
output series points should be indistinguishable. To
provide artificial distinction, some plots were made
with the input time series points denoted by the letter
I and the offset and normalized output time series
points denoted by the letter 0. Figure 11 is an example
of one of these plots.
With the calculated phase delay verified, a time
domain measure of the Modem system fidelity could be
calculated. Let Signal be defined as the squared value
of the input cosine time series at a given time index.
Let Noise be defined as the square of the difference

Figure 11
Modem Simulation Phase Delay Verification

between the input time series value and the offset and
normalized output time series value at each time index.
Finally, let SNR be defined as the ratio of average
Signal to average Noise calculated over a representative
length of time series indexes and expressed on a decibel
scale. The Modem program adds up Signal values and
Noise values over a 1000 sample loop and calculates the
value of ten times the logarithm of their ratio. The
resulting SNR is written to the screen.
The final block of code in the demodulator simu-
lation program is a loop that calls the FFT subroutine
SPFFTC. The power spectrum of the decimated time series
is calculated from the complex gain returned by SPFFTC.
An example of a Modem output power spectrum is shown in
Figure 12.

Figure 12
Modem Simulation Output Spectrum Example

When the frequency modulating baseband signal is
a single cosine, two methods may be used to calculate
the integral portion of the continuous time FM equation
as shown in Appendix A. The integral of the single
cosine may either be mathematically evaluated, which
yields the sine of the argument, or approximated by a
sum of cosines. The FM1C0S.F0R simulation program
calculates each point in the time series using both
methods. Any disagreement between the two methods is
quantified by calculating the average squared difference
between the methods across the time series. Because the
first method is mathematically exact and the second is
an approximation, this average squared difference is
called the error. The error, expressed in decibels, is
written to the screen.
Table 2 displays the error values for several
combinations of baseband and maximum deviation frequen-
cies. In Table 2, the sample frequency and IF were
arbitrarily chosen to be 1MHz and 250kHz, respectively,
and the time series contained 8192 samples. The fre-

quencies were assumed to be large enough to have neglig-
ible influence on the results. The selected time series
length is derived below.
Baseband, Hz
100 500 lk 5k 10k 15k
10 -93 -91 -91 -91 -91 -91
Fd 100 -73 -71 -71 -71 -71 -71
Hz lk -53 -51 -51 -51 -51 -51
10k -33 -31 -31 -31 -31 -31
50k -19 -17 -17 -17 -17 -17
Table 2
Cosine Integration Method Error Values
It is clear from the data in Table 2 that the
modulation frequency has very little effect on the
difference in values calculated by the two cosine inte-
gration methods. The variation in error with Fd, how-
ever, appears to be quite predictable and linear. In
fact, every decade increase in Fd corresponds to a
decade decrease in the magnitude of the error. (Recall
that the plotted error values are based on the squared
difference.) The empirical equation corresponding to
the data is
K = Fd x (error magnitude),
where K is a constant equal to approximately 360,000.
In an attempt to decrease the error in the

summation approximation method of integral evaluation,
an additional statement was added to FM1C0S.F0R that
calculated the value of the integral using the trape-
zoidal rule of numerical integration. A second error
value, analogous to the error term described above, is
calculated for the trapezoidal rule method and written
to the screen. The program was rerun using the input
values shown in Table 2. Rather than improving on the
error performance of the previous approximation method,
the trapezoidal rule method resulted in errors that were
one to two decibels worse.
In order for the digital demodulator to provide
good performance, the sample rate must be chosen to
support the highest significant frequency component of
the input FM spectrum. Therefore, the demodulator
optimization phase began with an investigation of the
bandwidth of this spectrum. Three initial assumptions
were made: firsts the audio baseband extends from 50Hz
to 15kHz; second, the maximum deviation frequency, Fd,
is fixed at 75kHz to correspond to the commercial FM
radio broadcasting standard; and third, the expected IF
bandwidth of a single commercial FM station extends out
75kHz on either side of the center frequency.
With these assumptions in mind, and ignoring for
a moment the CD compatibility issue, an initial sample
rate (Fs) of 400kHz was chosen. The desire to center

the IF spectrum between OHz and (Fs/2)Hz led to a choice
of 100kHz for the IF. (Note: for ease of explanation,
the intermediate frequency, or IF, will hereafter be
referred to as the carrier frequency, or Fc.) Figure 13
is the FM spectrum resulting from use of the above para-
meters with a single baseband frequency (Fm) of 15kHz.
Figure 13 clearly shows that the FM bandwidth
corresponding to the chosen parameter values is well
over the expected 75kHz. A design based on these values
would suffer from distortion due to aliasing and the
exclusion of significant FM frequency bands. The search
for better performance led to the generation of more
graphs with higher sampling frequencies. The goal was
to find sample and carrier frequencies for which all
significant FM spectral lines were contained in the
middle 75 to 80 percent of the graph. In all cases, Fc
was chosen to be equal to 0.25FS to center the FM spec-
trum between OHz and (Fs/2)Hz.
Figure 14 appears to demonstrate that a sample
frequency of 800kHz satisfies the goal. However, the
graph in Figure 14 suffers from an elevated noise floor
which is an artifact caused by truncation error in the
time series input to the FFT procedure. The elevated
noise floor can cause masking of low level frequency
components that might otherwise be visible in the FFT.
This error can be minimized by careful selection

flHpiitmt vi. Frequency
HDD FREQ: 15000 NOD HAS: 1.0000 SETA: 5.00 HAH DEI) FREQ: 75000 SANPLE FREQ: 500000 IDEA SIMPLES
Figure 13
FM Spectrum for Fs = 400kHz, Fc = 100kHz
Fd = 75kHz, and F = 15kHz

- anpiitmt vs. Frequency
FM Spectrum for Fs = 800kHz, Fc = 200kHz
Fd = 75kHz, and Fm = 15kHz

of the sample, carrier, and baseband frequencies.
First, given a fixed number of steps in the time series,
Fs should be chosen to give an integer frequency step
size in the resulting FFT. Second, Fc needs to be an
integer multiple of this step size. Third, the time
series should span an integer number of cosine periods
for each baseband frequency in the FFT. Finally, the
lowest baseband frequency that can be accommodated is
given by the reciprocal of the length of time spanned by
the time series, which in turn is equal to the number of
samples in the series times the reciprocal of the sample
frequency. This can be written as
FMin = VN
These interrelated tradeoffs affected the final
choices for time series length, sample and carrier fre-
quencies, and minimum baseband frequency. As shown in
Figure 14, a sample frequency of about 800kHz appeared
to accommodate the FM spectrum resulting from the stated
parameters. The FFT procedure required that the number
of samples in the time series (N) be equal to an integer
power of two. For N = 8192, a sample frequency of
802,816Hz, a carrier frequency of 200,704Hz, and a
minimum baseband frequency of 98Hz meet all of the
requirements. The FM spectrum resulting from these
parameters is shown in Figure 15. Note that 12,544Hz is
an integer multiple of the minimum baseband frequency.

Audiitijdt gs. Freamncv

Figure 15 reveals two important points: first,
that the technique for eliminating FFT truncation error
works; and second, that the sample frequency is still
too low. A less serious problem is that the maximum
appropriate baseband frequency, which must be a multiple
of 98Hz and should be about 15kHz to accurately repre-
sent the highest frequencies in the audio band, is
limited to about 12.5kHz.
Based on the information in Figure 15, it was
decided to raise the sample frequency to about 1MHz.
Given the time series length of 8192 samples, the actual
sample frequency was chosen to be 999,424Hz, which
allows minimum and maximum baseband frequencies of 122Hz
and 15,616Hz, respectively. The graph presented in
Figure 16 shows that this configuration almost meets the
FM spectrum goal. However* note the power contained in
the frequency components outside of the original single-
sided target FM bandwidth of 75kHz.
It seems clear that commercial radio broadcas-
ters must process their audio signal to artificially
limit the resulting FM bandwidth. At this juncture, it
was decided to investigate the maximum potential fidel-
ity of the digital demodulator given the unmodified
output of the modulator simulation. The conscientious
simulation of the frequency content of a real world FM
station could be pursued, time permitting, after the

ftHHitude us.

factors affecting the overall Modem program performance
potential were well understood. With that in mind, and
since the sample rate had already been increased twice,
the maximum deviation frequency was lowered to 50kHz.
The resulting FM spectrum, shown in Figure 17, has all
of its significant spectral lines in the middle 75 to 80
percent of the graph.
Now that the simulated FM bandwidth was known,
the I and Q channel FIR filter frequency response could
be more closely tailored to the heeds of the demodu-
lator. Referring to Figure 17, if the center frequency
were shifted down to 0Hz, as it is at the output of the
I and Q channel multipliers, a low pass filter that was
flat out to (437,248 + 15616 249856)Hz, or about
203kHz, would pass all of the FM spectral components
visible in Figure 17 unattenuated. Figure 18 shows the
frequency response of the revised I and Q channel FIR
filter design using a Ghebyshev window function.
The ratio and arctangent processes are insensi-
tive to changes in sample and maximum deviation fre-
quencies, so the next simulation element to be examined
was the delay variable in the finite difference block.
In order to make an informed first guess at an
appropriate value for the delay variable, called DD for
Differentiator Delay, a better understanding of the
arctangent block output was required. To provide this

AHfriitufie . frtouencv

ftnplitn-w. Frtquincy
Figure 18
Revised I and Q Channel FIR Lowpass Filter
Frequency Response

understanding, additional simulation program and plot
routine code was written to generate time domain graphs
of the arctangent block output time series.
Figure 19 is an example of the arctangent block
output for Fs = lMz, Fc = 250kHz, Fd = 45kHz, and a
baseband frequency of 15000Hz. The maximum deviation
frequency was decreased to 45kHz to contain the arctan-
gent output to within +/-tt. This artifice resulted in a
continuous graph, which, while not required for correct
demodulator operation, eases the task of analysis.
Recall the constraint imposed on the finite
difference block by the method of correcting for angle
differences crossing the quadrant II/quadrant III boun-
dary. This constraint required that the difference
being calculated in the finite difference block have a
magnitude less than tt. Examination of steepest portion
of the curve shown in Figure 19 reveals that the mag-
nitude of the finite difference block output exceeds w
when DD exceeds eleven.
This limit on DD was borne out in the simu-
lations. The curve in Figure 20 is a time domain base-
band output of the Modem program that was generated with
the same frequency values used for Figure 19 and a DD
value of 12. The effects of the quadrant II/quadrant
III correction process can clearly be seen.

Figure 19
Arctangent Block Output Time Series

Figure 20
Modem Simulation Time Domain Output with
Differentiator Delay Value Too Large

To verify the threshold value of DD, the Modem
program was run with the frequencies used to generate
Figure 20 and a DD value of eleven. The resulting time
series are shown in Figure 21. As expected, the quad-
rant II/quadrant III correction procedure did not cause
any distortion in this case.
Note that the DD value of eleven results in a
slight misalignment of the input and output time series.
This is because the demodulator phase delay calculation,
= (Oj+Oj+Oj)/2 ,
which was developed in the previous chapter, requires
the sum of FIR filter and finite difference block orders
to be even if the input and output time series are to be
superimposed. In the above example, 01 = 119 and 02 =
0, so 03 needs to be odd for super imposition to work.
Since 03 is one less than the difference in indices of
the time series elements being operated on by the finite
difference block, or DD-1, DD itself should be even.
Because eleven is not even, the phase delay term is not
an integer, and the time series are not aligned.
Although the pre-decimation FIR low pass filter
has been included in the block diagram of the digital
demodulator, its presence has yet to be vindicated.
Figure 22 is an FFT of the undecimated output of the
double precision Modem program showing the presence of
noise components out of the audio band.

Figure 21
Modem Simulation Time Domain Output with
Differentiator Delay Value Acceptable

Y5.Fr Figure 22
Modem Simulation Spectrum without
Pre-decimation Filter

To reduce the out of band noise components, a
201 coefficient pre-decimation filter was inserted in
the Modem program. The frequency response of this
filter, which uses a Chebyshev window function, is shown
in Figure 23.
As expected, the pre-decimation filter removed
most of the out of band noise. Figure 24 shows the
undecimated Modem output spectrum with the filter in
At this point, because the double precision
Modem program was working fairly well, it was decided to
collect time domain SNR values. The method of calcu-
lation of these values was described in the previous
chapter. Table 3 is a collection of SNR results for
various combinations of baseband frequencies. The
following parameter values were used to generate the
data in Table 3:
Fs = 999,424Hz
Fc = 249,856Hz
Fd = 50kHz
DD = 2 samples
These values and most of the frequencies listed in Table
3 were chosen to minimize FFT truncation error as ex-
plained previously. The 10kHz and 11kHz were chosen to
reveal whether or not intermodulation distortion was
occurring in the demodulator.

fiHt I ;t'joe vt.TrMutricv
-JO -
-60 -
-SO -
-105 -
150 -)---------1 i"11----- | I - |------- | i i I i
t 50000 100000 150000 200000 250000 300000 350000 500000 350000 500000
Figure 23
Frequency Response of Pre-decimation FIR Filter
Chebyshev Window Function

. ..s. .. AHpTiWQi 05. FKMUt'iCy

-JO .
50 i
-60 -
-00 -
110- |i
_ irk i

0 J1232 62161 93636 121920 156160 107392 210621 219056
Figure 24
Modem Simulation Spectrum with
Pre-decimation Filter

Baseband Frequency, Hz SNR, dB
244 170
488 164
976 164
1952 162
3904 183
7808 160
15616 277
3904, 15616 43
244, 488, 976, 1952, 3904, 7808, 15616 51
10000 172
iiooo 183
10000, 11000 70
Table 3
Modem Simulation Initial Time Domain SNR Values
Double Precision
It is clear from the data in Table 3 that the
simulation provides excellent fidelity for a single
baseband frequency. Obviously, however, multiple base-
band frequencies are being distorted. Figure 25 is an
FFT of the Modem output time series for the seven har-
monic input case from Table 3. The output time series
has been decimated by sixteen to provide better reso-
lution of the baseband frequencies.
Figure 25 shows that there are no appreciable in
band spurious components that could be responsible for
the poor SNR result. Furthermore, FFTs run on the
undecimated Modem output time series corresponding to
the seven harmonic input case revealed no significant
out of band spurious frequencies. The Modem simulation
was not suffering from frequency distortion.

00 H
(toolituae vs. Fr'j*ricy

3300 7001 11712 15(15 1J520
27320 31232
Figure 25
Modem Simulation Spectrum for Seven
Baseband Harmonics

With frequency domain problems essentially ruled
out, attention could be diverted to the time domain.
Figure 26 is a graph of the amplitude of 500 consecutive
time series points versus an arbitrary sample number.
Both input and output time series points are displayed,
with the output series values normalized and offset in
time as described in the previous chapter.
The cause of the poor SNR is clearly visible as
amplitude errors in the output time series. What is not
clearly visible is the cause of the amplitude errors,
one possible source is passband ripple in the FIR low-
pass filters. The I and Q channel filters are required
for correct digital demodulator operation, but the pre-
decimation filter is not. To test the passband ripple
hypothesis, the SNR value for the seven baseband har-
monics case was recalculated with the pre-decimation
filter removed. The new value was 83dB, a 32dB improve-
ment. This suggests that filter passband ripple is at
least partially responsible for the poor time domain SNR
There are FIR design programs available that
produce low pass filters with maximally flat frequency
responses. The next step in the pursuit of better time
domain SNR values would be to use such a program to
design the I and Q channel filters. Rather than pursue
that sort of optimization, it was decided to use the

Figure 26
Modem Simulation Time Domain Output For
Seven Baseband Harmonics

remaining research time to add arbitrary phase angles to
the baseband cosines and noise to the FM time series.
Before those two enhancements were added, the
Modem source code was backed down to single precision.
The triple digit SNR values reported in Table 3 sug-
gested that the digital demodulator concept was valid.
It was time to begin tailoring the simulation to more
attainable quantization levels. While single precision
floating point calculations are a long way from sixteen
bit, fixed point arithmetic, it was a step in the right
Table 4 displays the same data as Table 3,
except all variables and intrinsic functions are in
single precision.
Baseband Frequency, Hz SNR, dB
244 109
488 112
976 114
1952 112
3904 112
7808 109
15616 105
3904, 15616 43
244, 488, 976, 1952, 3904, 7808, 15616 51
10000 112
11000 114
10000, 11000 70
Table 4
Modem Simulation Initial Time Domain SNR Values
Single Precision

Although the single frequency SNR values are
much lower than those obtained using double precision,
the multiple frequency values are essentially the same.
This occurs because the amplitude errors swamp out the
difference in calculation precision. Figure 26 repre-
sents the single precision time series as accurately as
it represents the double precision time series that
generated it. As before, the major cause of the ampli-
tude distortion was thought to be FIR filter passband
ripple. The theory was rechecked by removing the pre-
decimation filter from the single precision Modem
simulation. The resulting SNR of 82dB is an improvement
of 31dB over the Figure obtained with the pre-decimation
filter in place. Once again, this strongly implicates
FIR filter passband ripple.
All of the SNR values up to this point have been
based on input cosines with relative phase angles of 0 *.
To more closely simulate even the most simple musical
signals requires that the cosines have arbitrary phase
relative to each other. The creation of pseudorandom
phase values for the Modem simulation baseband cosines
was discussed in the previous chapter. The data pre-
sented in Table 5 show the results of incorporating
pseudorandom phase values with the multiple frequency
entries in Table 4. All variables and intrinsic func-
tions used for Table 5 were single precision.

Baseband Frequency, Hz SNR, dB
3904, 15616 43
244, 488, 976, 1952, 3904, 7808, 15616 51
10000, 11000 70
Table 5
Modem Simulation Initial Time Domain SNR Values
Random Phase, Single Precision
As Table 5 shows, the addition of pseudorandom
phase has little effect on the performance of the single
precision simulation as measured by the time domain SNR.
Figure 27 shows 500 consecutive points of the undeci-
mated output and normalized time series superimposed on
the input time series. The small amplitude errors that
give rise to the poor SNR value are clearly visible. As
in the two previous cases, removal of the pre-decimation
FIR filter improved the SNR value, in this case to 82dB.
Figure 28 depicts the FFT of the time series
shown in Figure 27 after it has been decimated by six-
teen. The raised noise floor compared to Figure 25 is
due to the single precision versus double precision
calculations. Other than the raised noise floor, there
are no in band spurious frequencies visible in this
case. The addition of the pseudorandom phase angles has
not degraded the simulation performance in the frequency
The final modification made to the Modem simu-
lation program was the addition of pseudorandom noise to

Figure 27
Modem Simulation Time Domain Output for Seven
Harmonics, Random Phase, Single Precision

flhd I rtuot iii Frtqutncy
Figure 28
Modem Simulation Spectrum for Seven Harmonics,
Random Phase, Single Precision

the FM time series. The peak level of the added noise
can be selected by the user as described in Chapter IV.
The input noise was defined to be the average squared
value of the pseudorandom noise. The input signal was
defined to be the average squared value of the FM time
series. The input SNR is the ratio of average signal
over average noise expressed on a decibel scale.
Figure 29 shows the FFT of the FM time series
corresponding to the set of seven baseband frequencies
used above with added pseudorandom noise of 0.01 peak
amplitude. The calculations were tried with both single
and double precision variables and functions with vir-
tually no difference in results. The input and output
SNRs were calculated to be 48dB and 46dB, respectively.
Figure 29 reveals the main problem with the
added noise: it's not bandlimited. As such, it does
not accurately reflect the signal that would be present
at the input of the digital demodulator. The actual
input signal is filtered before it reaches the demodu-
lator as shown in Figure 1. Thus, the noise power
should have rolled off appreciably at the frequency
extremes shown in Figure 29.
This rolloff could be created by inserting an
FIR bandpass filter ahead of the I and Q channel multi-
pliers. The presence of this filter would further
aggravate the time domain SNR values by virtue of its
passband ripple. For this reason and because it

150 i i i i t (| .
0 52050 120)20 1173)2 208055 312320 371700 037200 033712
Frequency in Hr
Figure 29
FM Time Series Spectrum for Seven Harmonics,
Random Phase, Random Noise

represented an element not present in the digital demo-
dulator block diagram, it was decided not to pursue the
passband filter option.
Even though it was not a realistic simulation,
the Modem program performance in the presence of the
added wideband noise was examined. Figure 30 is the FFT
of the output time series, decimated by sixteen, that
corresponds to the input signal shown in Figure 29. The
output noise floor is seen to lie somewhat below that of
the input, but not by much.
Table 6 summarizes the SNR data for several runs
of the single precision Modem program using the para-
meter values listed for Table 3. All frequencies had
random phase terms and the maximum random noise ampli-
tude was set to 0.01, which gave a calculated time
domain SNR value of 48dB. The data in Table 6 confirm
that the Modem simulation output SNR is approximately
equal to the input SNR in this case.
To see how the digital demodulator responded to
a noisier input signal, the simulation was run using all
of the parameter values listed for Table 6 with the
exception of the noise peak amplitude, which was
increased to 0.1. The simulation output data corre-
sponding to this input time domain SNR of 28dB are
presented in Table 7. Once again, the output SNR values
are essentially the same as the input values.

AhdI iFudt u$. Frtqutncy
Figure 30
Modem Simulation Spectrum for Seven Harmonics,
Random Phase, Single Precision, Random Noise

Baseband Frequency, Hz SNR, dB
244 48
488 49
976 49
1952 49
3904 49
7808 49
15616 49
3904, 15616 42
244, 488, 976, 1952, 3904, 7808, 15616 46
10000 49
11000 49
10000, 11000 48
Table 6
Modem Simulation Initial Time Domain SNR Values
Single Precision, 0.01 Peak Input Noise
Baseband Frequency, Hz SNR, dB
244 28
488 29
976 29
1952 29
3904 29
7808 29
15616 29
3904, 15616 28
244, 488, 976, 1952, 3904, 7808, 15616 28
10000 29
11000 29
10000, 11000 28
Table 7
Modem Simulation Initial Time Domain SNR Values
Single Precision, 0.1 Peak Input Noise

The data presented in the previous chapter
proves that the digital demodulator works quite well,
even with a noisy input signal. All of the reported
data was generated with somewhat arbitrary simulation
parameters for the FIR filters and the differentiator
delay. The demodulator performance would doubtlessly be
enhanced by an investigation into the optimum values for
these parameters.
As a case in point, the output SNR for the seven
frequency input case improved significantly when the
pre-decimation filter was deleted, indicating that
filter passband ripple may be degrading demodulator
performance. However, the output SNR for the 15616Hz
input case decreased when the pre-decimation filter was
removed because of out of band noise. A thorough anal-
ysis of FIR filters would most likely turn up a can-
didate that would provide satisfactory out of band
attenuation with minimal passband ripple. A similar
analysis might produce an I and Q channel filter that
provided enhanced demodulator response. Following the

incorporation of better filters, the differentiator
delay value could optimized over the audio band. It is
likely that these parameters are interdependent and that
the best overall configuration would have to be found by
an iterative procedure involving all three.
Finally, the relatively limited set of quanti-
zation values offered by digital hardware needs to be
considered. The effects of this limitation on the
performance of the digital demodulator remain to be
simulated. But, based on the results of this simulation
and the performance available from the typical 16-bit
quantization CD player, it is predicted that an FM radio
receiver using 16-bit quantization and providing an SNR
of 70dB to 80dB or better will be commercially available
within the next few years. By the end of the 1990's the
analog FM radio will be as rare as the LP record is now.

1. Stearns, Samuel D., and David, Ruth A. Signal
Processing Algorithms. New Jersey: Prentice-Hall,
2. Carlson, A. Bruce. Communication Systems: An Intro-
duction to Signals and Noise in Electrical Communi-
cation. 3rd ed. New York: McGraw-Hill, 1986.
3. IEEE Acoustics, Speech, and Signal Processing
Society; Digital Signal Processing Committee, ed.
Programs for Digital Signal Processing. New York:
IEEE Press, 1979.

Definitions of variables used in the derivation.
y(t) = carrier waveform
Ac = maximum carrier amplitude
8c(t) total instantaneous phase angle
f(t) = carrier frequency
8(t) = relative instantaneous phase angle
f,-(t) = instantaneous frequency deviation
fd = maximum frequency deviation constant
m(t) = message waveform
60 = phase angle at t = 0
1 = discrete time index variable
L = current value of 1
Ts = sampling interval
fm = frequency of single message tone

Y(t) = Re{Acexp(j6c(t))}
y (t) = Accos[9c(t) ]
0c(t) = 27rfct + 0 (t)
f,(t) = d_[0(t)]
27T dt
0(t) = 2?r
f,(t) = fjm(t)
0(t) = 27T *t fdm(r)d t -00
0(t) rt = 27rfd m(r)dr J-CO
0(t) rt = 27rfd m(r)dr + 0O Jo
= 0.
0(t) rt = 27rfd m(r)dr Jo
y(t) = Acos[27rfct + 27rfd m(T)dr]
Converting to discrete time,
y(LTs) Accos [ 27rfcLTs + 27rfd 2 m(lTs)

For single tone modulation, where m(t) = cos(27rfmt),
y(t) Accos[27rfct + 2n"fd cos (2?rfn]TdT) ]
y (t) = AcCOS[27Tfct + fj sin(27Tfmt) ]
Converting to discrete time,
y(LTs) = Accos[27rfcLTs + fdsin(27rfmlTs) ]
y(LTs) Accos[27TfcLTs + 27rfdTs 2 cos (27rfmlT8) ]

Program FM1COSDP
C ** Program to frequency modulate a single baseband **
C ** cosine. Sample, modulating, and maximum ****
C ** deviation frequencies entered via keyboard. *****
C ** Output is 512 point FFT written to ASCII file. ***
Double Precision TwoPi, OmegaMt, Ts, OmegaCt, DfTs,
* Fs, Beta, XI, X2, X3, Error1, Error2, Mag,
* AmpMax, Fc, FCos, FDelta, Suml, Sum2
Complex*16 X(0:8200)
Dimension Mag(0:512)
Open (unit=l, file='C:\RMF0R\FM1DB.DAT)
Write(*,*) ('Enter sample frequency: ')
Read(*,*) Fs
Write(*,*) ('Enter modulating frequency: ')
Read(*,*) FCos
Write(*,*) ('Enter Max deviation frequency: ')
Read(*,*) FDelta
C ** Establish constants, initialize variables *********
N = 1024
Am = 1.D0
Suml = 0.D0
Sum2 = 0.D0
Errorl = 0.D0
Error2 = 0.D0
Fc = Fs/4.DO
Ts = l.DO/Fs
TwoPi = 8.D0*Atan(l.D0)
OmegaMt = TwoPi*FCos*Ts
OmegaCt = TwoPi*Fc*Ts
DfTs = TwoPi*FDelta*Ts
Beta = FDelta/FCos
C ** Calculate FM'ed single cosine time series *********
C ** using three methods. Find errors in two *********
C ** summation approximation methods relative to ******
C ** exact method. Write errors to screen in dB. *****