GaussianWaves

Sampling Theorem - Bandpass or Intermediate or Under Sampling

2011-07-28T11:40:00.001+01:00

Intermediate Sampling or Under-Sampling

If the signal is a bandpass signal with bandwidth Fb (Bandwidth is simply the difference between the lowest and the highest frequency present in the signal).

In order for a faithful reproduction and reconstruction of a bandpass analog signal with bandwidth - Fb, the signal should be sampled at a Sampling frequency (Fs) that is greater than or equal to twice the maximum bandwidth of the signal.

Consider a bandpass signal extending from 150Hz to 200Hz. The bandwidth of this signal is Fb=200Hz-150Hz=50Hz. In order to faithfully represent the above signal in the digital domain the sampling frequency must be Fs>=100Hz (2*Fb). Note that the sampling frequency 100Hz is far below the maximum content of the signal (which is 200Hz). That is why the bandpass sampling is also called "under-sampling". As long as the sampling frequency is greater than or equal to twice the bandwidth of the signal, the reconstruction back to analog domain will be error free.

Going back to the aliasing zone figure, if the signal of interest is in the zone other than zone 1, it is called a bandpass signal and the sampling operation is called "Intermediate Sampling" or "Harmonic Sampling" or "Under Sampling" or "Bandpass Sampling".

Note that zone 1 is a mirror image of zone 2 (with frequency reversal). Similarly zone 3 is a mirror image of zone 4 etc.., Also, any signal in zone 1 will be reflected in zone 2 with frequency reversal which inturn will be copied in zone 3 and so on.

Lets say the signal of interest lies in zone 2. This will be copied in all the other zones. Zone 1 also contains the sampled signal with frequency reversal which can be correct by reversing the order of FFT in digital domain.

No matter in which zone the signal of interest lies, zone 1 always contains the signal after sampling operation is performed. If the signal of interest lies in any of the even zones, zone 1 contains the sampled signal with frequency reversal. If the signal of interest lies in any of the odd zones, zone 1 contains the sampled signal without frequency reversal.

Example:

Consider an AM signal centered at carrier frequency 1MHz, with two components offset by 10KHz - 0.99 MHz and 1.01 MHz. So the AM signal contains three frequency components at 0.99 MHz, 1 MHz and 1.01 MHz.

Our desire is to sample the AM signal. As with the usual sampling theorem (baseband), we know that if we sample the signal at twice the maximum frequency i.e Fs>=2*1.01MHz=2.02 MHz there should be no problem in representing the analog signal in digital domain.

By the bandpass sampling theorem, we do not need to use a sampler running at Fs>=2.02 MHz. Faster sampler implies more cost. By applying the bandpass sampling theorem, we can use a slower sampler and reduce the cost of the system. The bandwidth of the signal is 1.01MHz-0.99 MHz = 20 KHz. So, just sampling at Fs>=2*20KHz=40KHz will convert the signal to digital domain properly and we can also avoid using an expensive high rate sampler (if Fs >= 2.02 MHz used according to baseband sampling theorem).

Lets set the sampling frequency to be Fs=120KHz (which is 3 times higher than the minimum required sampling rate of 40KHz or oversampling rate =3).

Now we can easily find the position of the spectral components in the sampled output by using the aliasing zone figure as given above. Since Fs=120 KHz, Fs/2 will be 60KHz. So the zone 1 will be from 0 to 60 KHz, zone 2 -> 60-120KHz and so on. The three spectral components at 0.99MHz, 1MHz and 1.01 MHz will fall at zone 17 ( how ? 0.99 MHz/60 KHz = 16.5 , 1MHz/60KHz = 16.67 and 1.01MHz/60KHz = 16.83 , all figures approximating to 17). By the aliasing zone figure, zone 16 contains a copy of zone 17, zone 15 contains a copy of zone 16, zone 14 contains a copy of zone 15 and so on... Finally zone 1 contains the copy of zone 2 (Frequency reversal also exist at even zones). In effect, zone 1 contains a copy of zone 17. Since the original spectral components are at zone 17, which is an odd zone, zone 1 contains the copy of spectral components at zone 17 without frequency reversal.

Since there is no frequency reversal, in zone 1 the three components will be at 30KHz, 40KHz and 50KHz (You can easily figure this out ).

This operation has downconverted our signal from zone 17 to zone 1 without distorting the signal components. The downconverted signal can be further processed by using a filter to select the baseband downconverted components. Following figure illustrates the concept of bandpass sampling.

Consider the same AM signal with three components at 0.99MHz, 1MHz and 1.01 MHz. Now we also have an "unwanted" fourth component at 1.2 MHz along with the incoming signal. If we sample the signal at 120 KHz, it will cause aliasing (because the bandwidth of the entire signal is 1.2-0.99 = 0.21 MHz = 210 KHz and the sampling frequency of 120KHz is below twice the bandwidth). In order to avoid anti-aliasing and to discard the unwanted component at 1.2 MHz, an anti-aliasing bandpass filter has to be used to select those desired component before performing the sampling operation at 120KHz. This is also called "pre-filtering". The following figure illustrates this concept.

Recommended Books:

Sampling Theorem - Baseband Sampling

2011-07-21T10:04:00.000+01:00

"Nyquist-Shannon Sampling Theorem" is the fundamental base over which all the digital processing techniques are built.

Processing a signal in digital domain gives several advantages (like immunity to temperature drift, accuracy, predictability, ease of design, ease of implementation etc..,) over analog domain processing.

Analog to Digital conversion:

In analog domain, the signal that is of concern is continuous in both time and amplitude. The process of discretization of the analog signal in both time domain and amplitude levels yields the equivalent digital signal.

The conversion of analog to digital domain is a three step process
1) Discretization in time - Sampling
2) Discretization of amplitude levels - Quantization
3) Converting the discrete samples to digital samples - Coding/Encoding

The sampling operation samples ("chops") the incoming signal at regular interval called "Sampling Rate" (denoted by Ts). Sampling Rate is determined by Sampling Frequency (denoted by Fs) as

Lets consider the following logical questions:
* Given a real world signal, how do we select the sampling rate in order to faithfully represent the signal in digital domain ?
* Is there any criteria for selecting the sampling rate ?
* Will there be any deviation if the signal is converted back to analog domain ?

Answer : Consult the "Nyquist-Shannon Sampling Theorem" to select the sampling rate or sampling frequency.

Nyquist-Shannon Sampling Theorem:
The following sampling theorem is the exact reproduction of text from Shannon's classic paper[1],

"If a function f(t) contains no frequencies higher than W cps, it is completely determined by giving its ordinates at a series of points spaced 1/2W seconds apart."

Sampling Theorem mainly falls into two categories :

1) Baseband Sampling - Applied for signals in the baseband (useful frequency components extending from 0Hz to some Fm Hz)

2) Bandpass Sampling - Applied for signals whose frequency components extent from some F1 Hz to F2Hz (where F2>F1)

In simple terms, the Nyquist Shannon Sampling Theorem for Baseband can be explained as follows

Baseband Sampling

If the signal is confined to a maximum frequency of Fm Hz, in other words, the signal is a baseband signal (extending from 0 Hz to maximum Fm Hz).

In order for a faithful reproduction and reconstruction of an analog signal that is confined to a maximum frequency Fm, the signal should be sampled at a Sampling frequency (Fs) that is greater than or equal to twice the maximum frequency of the signal.

Consider a 10Hz sine wave in analog domain. The maximum frequency present in this signal is Fm=10 Hz (obviously no doubt about it !!!). Now, to satisfy the sampling theorem that is stated above and to have a faithful representation of the signal in digital domain, the sampling frequency can be chosen as Fs >=20Hz. That is, we are free to choose any number above 20 Hz. Higher the sampling frequency higher is the accuracy of representation of the signal. Higher sampling frequency also implies more samples, which implies more storage space or more memory requirements.

In time domain, the process of sampling can be viewed as multiplying the signal with a series of pulses ("pulse train) at regular interval - Ts. In frequency domain, the output of the sampling process gives the following components - Fm (original frequency content of the signal), Fs±Fm,2Fs±Fm,3Fs±Fm,4Fs±Fm and so on and on...

Now the sampled signal contains lots of unwanted frequency components (Fs±Fm,2Fs±Fm,...). If we want to convert the sampled signal back to analog domain, all we need to do is to filter out those unwanted frequency components by using a "reconstruction" filter (In this case it is a low pass filter) that is designed to select only those frequency components that are upto Fm Hz. The above process mentions only the sampling part which samples the incoming analog signal at regular intervals. Actually a quantizer will follow the sampler which will discretize ("quantize") amplitude levels of the sampled signal. The quantized amplitude levels are sent to an encoder that converts the discrete amplitude levels to binary representation (binary data). So when converting the binary data back to analog domain, we need a Digital to Analog Converter (DAC) that converts the binary data to analog signal. Now the converted signal after the DAC contains the same unwanted frequencies as well as the wanted component. Thus a reconstruction filter with proper cut-off frequency has to placed after the DAC to filter out only the wanted components.

Aliasing and Anti-aliasing:

Consider a signal with two frequency components f1=10Hz - which is our desired signal and f2=20Hz - which is a noise.Lets say we sample the signal at 30Hz. The first frequency component f1=10Hz will generate following frequency components at the output of the multiplier (sampler) - 10Hz,20Hz,30Hz,50Hz,70Hz and so on. The second frequency component f2=20Hz will generate the following frequency components at the output of the multiplier - 20Hz,10Hz,50Hz,40Hz,80Hz and so on...

Note the 10Hz component that is generated by f2=20Hz. This 10Hz component (which is a manifestation of noisy component f2=20Hz) will interfere with our original f1=10Hz component and are indistinguishable.This 10Hz component is called "alias" of the original component f2=20Hz (noise). Similarly the 20Hz component generated by f1=10Hz component is an "alias" of f1=10Hz component. This 20Hz alias of f1=10Hz will interfere with our original component f2=20Hz and are indistinguishable. We do not need to care about the interference that occurs at 20Hz since it is a noise and any way it has to be eliminated. But we need to do something about the aliasing component generated by the f2=20Hz. Since this is a noise component, the aliasing component generated by this noise will interfere with our original f1=10Hz component and will corrupt it.

Aliasing depends on the sampling frequency and its relationship with the frequency components. If we sample a signal at Fs, all the frequency components from Fs/2 to Fs will be alias of frequency components from 0 to Fs/2 and vice versa. This frequency - Fs/2 is called "Folding frequency" since the frequency components from Fs/2 to Fs folds back itself and interferes with the components from 0Hz to Fs/2 Hz and vice versa. Actually the aliasing zones occur on the either sides of 0.5Fs, 1.5Fs, 2.5Fs,3.5Fs etc... All these frequencies are also called "Folding Frequencies" that causes frequency reversal. Similarly aliasing also occurs on either side of Fs,2Fs,3Fs,4Fs... without frequency reversals. The following figure illustrates the concept of aliasing zones.

In the above figure, zone 2 is just a mirror image of zone 1 with frequency reversal. Similarly zone 2 will create aliases in zone 3 (without frequency reversal), zone 3 creates mirror image in zone 4 with frequency reversal and so on...

In the example above, the folding frequency was at Fs/2=15Hz, so all the components from 15Hz to 30Hz will be the alias of the components from 0Hz to 15Hz. Once the aliasing components enter our band of interest, it is impossible to distinguish between original components and aliased components and as a result, the original content of the signal will be lost. In order to prevent aliasing, it is necessary to remove those frequencies that are above Fs/2 before sampling the signal. This is achieved by using an "anti-aliasing" filter that precedes the analog to digital converter.

An anti-aliasing filter is designed to restrict all the frequencies above the folding frequency Fs/2 and therefore avoids aliasing that may occur at the output of the multiplier otherwise.

Thus, a complete design of analog to digital conversion contains an anti-aliasing filter preceding the ADC and the complete design of digital to analog conversion contains a reconstruction filter succeeding the DAC.

Note: Remember that both the anti-aliasing and reconstruction filters are analog filters since they operate on analog signal. So it is imperative that the sampling rate has to be chosen carefully to relax the requirements for the anti-aliasing and reconstruction filters.

Effects of Sampling Rate:

Consider a sinusoidal signal of frequency Fm=2MHz. Lets say that we sample the signal at Fs=8MHz (Fs>=2*Fm). The factor Fs/Fm is called "over-sampling factor". In this case we are over-sampling the signal by a factor of Fm=8MHz/2MHz = 4. Now the folding frequency will be at Fs/2 = 4MHz. Our anti-aliasing filter has to be designed to strictly cut off all the frequencies above 4MHz to prevent aliasing.

In practice, ideal brick wall response for filters is not possible. Any filter will have a transition band between pass-band and stop-band. Sharper/faster roll off transition band (or narrow transition band) filters are always desired. But such filters are always of high orders. Since both the anti-aliasing and reconstruction filters are analog filters, high order filters that provide faster roll-off transition bands are expensive (Cost increases proportionally with filter order). The system also gets bulkier with increase in filter order.Therefore, to build a relatively cheaper system, the filter requirement in-terms of width of the transition band has to be relaxed. This can be done by increasing the sampling rate or equivalently the over-sampling factor. When the sampling rate (Fs) is increased, the distance between the maximum frequency content Fm and Fs/2 will increase. This increase in the gap between the maximum frequency content of the signal and Fs/2 will ease requirements on the transition bands of the anti-aliasing analog filter.

Following figure illustrates this concept. If we use a sampling frequency of Fs=8MHz (over-sampling factor = 4), the transition band is narrower and it calls for a higher order anti-aliasing filter (which will be very expensive and bulkier). If we increase the sampling frequency to Fs=32MHz (over-sampling factor = 32MHz/2MHz=16), the distance between the desired component and Fs/2 has greatly increased that it facilitates a relatively inexpensive anti-aliasing filter with a wider transition band.

Thus, increasing the sampling rate of the ADC facilitates simpler lower order anti-aliasing filter as well as reconstruction filter. However, increase in the sampling rate calls for a faster sampler which makes ADC expensive. It is necessary to compromise and to strike balance between the sampling rate and the requirement of the anti-aliasing/reconstruction filter.

Application : Up-conversion

In the above examples, the reconstruction filter was conceived as a low pass filter that is designed to pass only the baseband frequency components after reconstruction. Remember that any frequency component present in zone 1 will be reflected in all the zones (with frequency reversals in even zones and without frequency reversals in odd zones). So, if we design the reconstruction filter to be a bandpass filter selection the reflected frequencies in any of the zones expect zone 1, then we achieve up-conversion.

In any communication system, the digital signal that comes out of a digital signal processor cannot be transmitted across as such. The processed signal (which is in the digital domain) has to be converted to analog signal and the analog signal has to be translated to appropriate frequency of operation that fits the medium of transmission. For example, in an RF system, a baseband signal is converted to higher frequency (up-conversion) using a multiplier and oscillator and then the high frequency signal is transmitted across the medium. If we have a band-pass reconstruction filter at the output of the DAC, we can directly achieve up-conversion (which saves us from using a multiplier and oscillator). The following figure illustrates this concept.

Reference:

[1] Claude E. Shannon, "Communication in the presence of noise" ,Proceedings of the IRE, Vol 37, no.1,pp.10-21,Jan 1949

Recommended Books:

Simulation of OFDM system in Matlab - BER Vs Eb/N0 for OFDM in AWGN channel

2011-07-08T11:16:00.005+01:00

To simulate an OFDM system, following design parameters are essential.
Lets consider the OFDM system parameters as defined in IEEE 802.11^[1] specifications

Given parameters in the spec:
N=64; %FFT size or total number of subcarriers (used + unused) 64
Nsd = 48;    %Number of data subcarriers 48
Nsp = 4 ;    %Number of pilot subcarriers 4
ofdmBW = 20 * 10^6 ; % OFDM bandwidth

Derived Parameters:
deltaF = ofdmBW/N;      % Bandwidth for each subcarrier - include all used and unused subcarries
Tfft = 1/deltaF;               % IFFT or FFT period = 3.2us
Tgi = Tfft/4;                    % Guard interval duration - duration of cyclic prefix - 1/4th portion of OFDM symbols
Tsignal = Tgi+Tfft;    % Total duration of BPSK-OFDM symbol = Guard time + FFT period
Ncp = N*Tgi/Tfft; %Number of symbols allocated to cyclic prefix
Nst = Nsd + Nsp;           % Number of total used subcarriers
nBitsPerSym=Nst; %For BPSK the number of Bits per Symbol is same as number of subcarriers

We will use the following model for simulating the OFDM system.

The channel in this case is modeled as a simple AWGN channel. Since the channel is considered to be an AWGN channel, there is no need for the frequency domain equalizer in the OFDM receiver ( Frequency domain equalizer will be helpful only if the channel introduces multipath fading). Since our channel is an AWGN channel, the frequency domain equalizer block in the above diagram can be removed.

Calculating Es/N0 or Eb/N0 for OFDM system:
In order to do a Monte Carlo simulation of an OFDM system, required amount of channel noise has to be generated that is representative of required Eb/N0. In Matlab it is easier to generate a Gaussian noise with zero mean and unit variance. The generated zero-mean-unit-variance noise has to be scaled accordingly to represent the required Eb/N0 or Es/N0. If we have Es/N0, the required noise can be generated from zero-mean-unit-variance-noise by,

Since the OFDM system transmits and received the data in symbols, it is appropriate/easier to generate required noise based on Es/N0 instead of Eb/N0 (as given above). But we are interested in plotting BER against Eb/N0. Ok!!! Then... how do we convert given Eb/N0 to Es/N0 for an OFDM system ?

Normally for a simple BPSK system, bit energy and symbol energy are same. So Eb/N0 and Es/N0 are same for a BPSK system. But for a OFDM BPSK system, they are not the same. This is because, each OFDM symbol contains additional overhead in both time domain and frequency domain. In the time domain, the cyclic prefix is an additional overhead added to each OFDM symbol that is being transmitted. In the frequency domain, not all the subcarriers are utilized for transmitted the actual data bits, rather a few subcarriers are unused and are reserved as guard bands.

Effect of Cyclic Prefix on Es/N0:
The following diagram illustrates the concept of cyclic prefix. Each OFDM symbol contains both useful data and overhead (in the form of cyclic prefix). The bit energy represents the energy contained in the useful bits. In this case, the bit energy is spread over N bits (where N is the FFT size). On top of the useful data, additional Ncp bits are added as cyclic prefix, which forms the overhead. So if the entire OFDM symbol is considered, the symbol energy is spread across N+Ncp bits.

This relationship is given as

which translates to,

Effect of unused subcarriers on Es/N0:

As mentioned earlier, not all the subcarriers are used for transmission. Out of the total N subcarriers, only Nst carriers are used for OFDM symbols transmission (this includes both data and pilot subcarriers). Again, in the frequency domain, the useful bit energy is spread across Nst subcarriers, whereas the symbol energy is spread across N subcarriers. This gives us another relationship between Es and Eb as given below

which translates to,

From (1) and (2), the overall effect of both cyclic prefix and unused subcarriers on Es/N0 is given by

which, when converted to dB yields the following relationship.

Since Ncp cylic prefix are added to the OFDM symbol, the output signal from the parallel to serial converter has to be boosted to compensate for the wastage of energy associated with the addition of cyclic prefix.
To properly generate the required SNR in Matlab, the signal term at the output of the parallel to serial converter has to be scaled as follows

The received signal is represented as (for the given Eb/N0)
y = boosted OFDM signal + required noise

Arrangement of subcarriers:

The IEEE 802.11 specification, specifies how to arrange the given subcarriers. The 52 used subcarriers (data + pilot) are assigned numbers from -26,-25,...-2,,-1 and 1,2,...,25,26. The following figure illustrates the scheme of assigning these subcarriers to the IFFT inputs.

Matlab Script: click here to download

Simulated Result:
From the simulated result, it can be ascertained that the OFDM- BPSK modulation has no advantage over a normal BPSK system in AWGN. OFDM proves to be effective in multipath environments (simulation of OFDM-BPSK system in a Rayleigh multipath environment will follow in future posts).

References
[1] IEEE 802.11 specification - "Orthogonal frequency division multiplexing (OFDM) PHY specification for the 5 GHz band" - chapter 17

More on OFDM:
Introduction to OFDM - part 1
Introduction to OFDM - part 2 - Basic OFDM transmitter
Introduction to OFDM - part 3 - OFDM System
Introduction to OFDM - part 4 - Cyclic Prefix
Simulation of OFDM System - Eb/N0 Vs BER for OFDM-BPSK over AWGN Channel

Recommended Books:

Introduction to OFDM - orthogonal Frequency division multiplexing - part 4 - Cyclic Prefix

2011-07-06T08:06:00.009+01:00

Consider a non-ideal channel h(t)≠δ(t), that causes delay dispersion. Delay dispersion manifests itself as Inter Symbol Interference (ISI)on each subcarrier channel due to pulse overlapping. It will also cause ICC (Inter Carrier Interference ) due to the non-orthogonality of the received signal. Adding cyclic prefix to each OFDM symbol mitigates the problems of ISI and ICC by removing them altogether.

Lets say, without cyclic prefix we transmit the following N values (N=Nfft=length of FFT/IFFT) for a single OFDM symbol.

Lets consider a cyclic prefix of length Ncp, ( where Ncp<N ), is formed by copying the last Ncp values from the above vector of X and adding those Ncp values to the front part of the same X vector.

With a cyclic prefix length Ncp, ( where Ncp<N ), the following values constitute a single OFDM symbol

If T is the duration of the an OFDM symbol in secs, due to the addition of cyclic prefix of length Ncp, the total duration of an OFDM symbol becomes T+Tcp, where Tcp=Ncp*T/N. Therefore, the number of samples allocated for cyclic prefix can be calculated as Ncp=Tcp*N/T, where N is the FFT/IFFT length, T is the IFFT/FFT period and Tcp is the duration of cyclic prefix.

The key ideas behind adding cyclic prefix :
1) Convert linear convolution in to circular convolution which eases the process of detecting the received signal by using a simple single tap equalizer
If you wish to know how the addition of cyclic prefix converts linear convolution to circular convolution, visit this link
2) Help combat ISI and ICC.

When a cyclic prefix of length Ncp is added to the OFDM symbol, the output of the channel (r) is given by circular convolution of channel impulse response (h) and the OFDM symbols with cyclic prefix (x).

As we know, for the discrete signals, circular convolution in the time domain translates to multiplication in the frequency domain.

Thus the above equation translates to

in the frequency domain.

At the receiver, R is the received signal (in Frequency domain) and our goal is to estimate the transmitted signal (X) from the received signal R. From the above equation, the problem of detecting the transmitted signal at the receiver side translate to a simple equalization equation as follows

After the FFT performed at the receiver side (i.e. after the FFT block in the receiver side), a single tap equalizer (which essentially implements the above equation) is used to estimate the transmitted OFDM symbol. It also corrects the phase and equalizes the amplitude.

A basic OFDM architecture with Cyclic Prefix is given below. (In the following diagram, symbols represented by small case letters are assumed to be in time domain, whereas the symbols represented by uppercase letters are assumed to be in frequency domain)

The IEEE specs specify the length of the Cyclic prefix in terms of its duration.
Lets see how to convert the specified duration (Tcp) in to actual number of samples assigned for cyclic prefix (Ncp).

Lets see an example of deriving Ncp from IEEE 802.11 spec ^[1]

Given parameters in the spec:
N=64; %FFT size or total number of subcarriers (used + unused) 64
Nsd = 48;    %Number of data subcarriers 48
Nsp = 4 ;    %Number of pilot subcarriers 4
ofdmBW = 20 * 10^6 ; % OFDM bandwidth

Derived Parameters:
deltaF = ofdmBW/N;      % Bandwidth for each subcarrier - include all used and unused subcarries
Tfft = 1/deltaF;               % IFFT or FFT period = 3.2us
Tgi = Tfft/4;                    % Guard interval duration - duration of cyclic prefix - 1/4th portion of OFDM symbols
Tsignal = Tgi+Tfft;    % Total duration of BPSK-OFDM symbol = Guard time + FFT period
Ncp = N*Tgi/Tfft; %Number of symbols allocated to cyclic prefix
Nst = Nsd + Nsp;           % Number of total used subcarriers
nBitsPerSym=Nst; %For BPSK the number of Bits per Symbol is same as number of subcarriers

References:
[1] IEEE 802.11 specification - "Orthogonal frequency division multiplexing (OFDM) PHY specification for the 5 GHz band" - chapter 17

More on OFDM:
Introduction to OFDM - part 1
Introduction to OFDM - part 2 - Basic OFDM transmitter
Introduction to OFDM - part 3 - OFDM System
Introduction to OFDM - part 4 - Cyclic Prefix
Simulation of OFDM System - Eb/N0 Vs BER for OFDM-BPSK over AWGN Channel

Introduction to OFDM - orthogonal Frequency division multiplexing - part 3

2011-06-27T11:03:00.003+01:00

In the previous post, the architecture of an OFDM transmitter was described using sinusoidal components. Generally, an OFDM signal can be represented as

s(t) = symbols mapped to chosen constellation (BPSK/QPSK/QAM etc..,)
Fn = orthogonal frequency

This equation can be thought of as an IFFT process ( Inverse Fast Fourier Transform). The Fourier transform breaks a signal into different frequency bins by multiplying the signal with a series of sinusoids. This essentially translates the signal from time domain to frequency domain. But, we always view IFFT as a conversion process from frequency domain to time domain.

FFT is represented by

where as its dual , IFFT is given by

The equation for FFT and IFFT differ by the co-efficients they take and the minus sign. Both equation does the same thing. They multiply the incoming signal with a series of sinusoids and separates them into bins.
In fact, FFT and IFFT are dual and behaves in a similar way.IFFT and FFT blocks are interchangeable.

Since the OFDM signal ( c(t) in the equation above ) is in time domain, IFFT is the appropriate choice to use in the transmitter, which can be thought of as converting frequency domain samples to time domain samples. Well, you might ask : s(t) is not in frequency domain and they are already in time domain; so whats the need to convert it into time domain again ? The answer is IFFT/FFT equation comes handy in implementing the conversion process and we can eliminate the individual sinusoidal multipliers required in the transmitter/receiver side. The following figure illustrates, how the use of IFFT in the transmitter eliminates the need for separate sinusoidal converters. Always remember that IFFT and FFT blocks in the transmitter are interchangeable as long as their duals are used in receiver.

The entire architecture of a basic OFDM system with both transmitter and receiver will look like this

An OFDM system is defined by IFFT/FFT length - N ,the underlying modulation technique ( BPSK/QPSK/QAM), supported data rate, etc..,. The FFT/IFFT length N defines the number of total subcarriers present in the OFDM system. For example, an OFDM system with N=64 provides 64 subcarriers. In reality, not all the subcarriers are utilized for data transmission. Some subcarriers are reserved for pilot carriers (used for channel estimation/equalization and to combat magnitude and phase errors in the receiver) and some are left unused to act as guard band. OFDM system do not transmit any data on the subcarriers that are near the two ends of the transmission band ( Not necessarily at the ends of the bands, implementation may differ). These subcarriers are collectively called guard band. The reservation of subcarriers to guard bands helps to reduce the out of band radiation and thus eases the requirements on transmitter front-end filters.The subcarriers in the guard band are also called Null subcarriers or virtual subcarriers.

For example IEEE 802.11 standard^[1] specifies the following parameters for its OFDM physical layer.

FFT/IFFT size = 64 ( implies 64 subcarriers in total = used + unused = Nfft)
Number of data subcarriers = 48 ( Nd)
Number of pilot subcarriers = 4 (Np)

Derived parameters from the above specification.
Number of total USED subcarriers = 52 ( Nu = Nd+ Np )
Number of total UNUSED subcarriers = 12 ( Nun = Nfft - Nu ).

According to the spec, these 52 used subcarriers are distributed in the following way. The 52 used subcarriers are named as 1,2,3,...,26 and -1,-2,-3,...,-26. The used subcarriers 1 to 26 are mapped to 1 to 26 of IFFT inputs and the subcarriers -1,-2,..,-26 are mapped to the IFFT inputs 38 to 63. The remaining IFFT inputs 27 to 37 and the input 0 (dc input) are set to 0 .In this manner the 12 null subcarriers are mapped to IFFT inputs.

[1] IEEE 802.11 specification - "Orthogonal frequency division multiplexing (OFDM) PHY specification for the 5 GHz band" - chapter 17

Introduction to OFDM - orthogonal Frequency division multiplexing - part 2

2011-06-08T11:09:00.006+01:00

Consider that we want to send the following data bits using OFDM : D = {d0,d1,d2,...). The first thing that should be considered in designing the OFDM transmitter is the number of subcarriers required to send the given data. As a generic case, lets assume that we have N subcarriers. Each subcarriers are centered at frequencies that are orthogonal to each other (usually multiples of frequencies).

The second design parameter could be the modulation format that we wish to use. An OFDM signal can be constructed using anyone of the following digital modulation techniques namely BPSK, QPSK, QAM etc..,
The data (D) has to be first converted from serial stream into parallel stream depending on the number of sub-carriers (N). Since we assumed that there are N subcarriers allowed for the OFDM transmission, we name the subcarriers from 0 to N-1. Now, the Serial to Parallel converter takes the serial stream of input bits and outputs N parallel streams (indexed from 0 to N-1). These parallel streams are individually converted into the required digital modulation format (BPSK, QPSK, QAM etc..,). Lets call this output S₀,S₁,..S_N. The conversion of parallel data (D) into the digitally modulated data (S) is usually achieved by a constellation mapper, which is essentially a look up table (LUT). Once the data bits are converted to required modulation format, they need to be superimposed on the required orthogonal subcarriers for transmission. This is achieved by a series of N parallel sinusoidal oscillators tuned to N orthogonal frequencies (f₀,f₁,...f_N-1). Finally, the resultant output from the N parallel arms are summed up together to produce the OFDM signal.

The following figure illustrates the basic concept of OFDM transmission (note: In order to give a simple explanation to illustrate the underlying concept,the usual IFFT/FFT blocks that are used in actual OFDM system, are not used in the block diagram)

OFDM Transmitter: (click to enlarge)

Example: (click images to enlarge)
The first example illustrates the concept of OFDM transmission with BPSK modulation as its underlying modulation format. The second example illustrates the OFDM transmission with pi/4 shifted QPSK modulation. Here 5 orthogonal subcarriers are assumed for the OFDM transmission.

Introduction to OFDM - orthogonal Frequency division multiplexing

2011-05-24T08:00:00.004+01:00

OFDM
In modulations, information is mapped on to changes in frequency, phase or amplitude (or a combination of them) of a carrier signal. Multiplexing deals with allocation/accommodation of users in a given bandwidth (i.e. it deals with allocation of available resource).

OFDM is a combination of modulation and multiplexing. In this technique, the given resource (bandwidth) is shared among individual modulated data sources. Normal modulation techniques (like AM, PM, FM, BPSK, QPSK, etc.., ) are single carrier modulation techniques, in which the incoming information is modulated over a single carrier. OFDM is a multicarrier modulation technique, which employs several carriers, within the allocated bandwidth, to convey the information from source to destination. Each carrier may employ one of the several available digital modulation techniques (BPSK, QPSK, QAM etc..,).

Why OFDM
OFDM is very effective for communication over channels with frequency selective fading ( different frequency components of the signal experience different fading). It is very difficult to handle frequency selective fading in the receiver , in which case, the design of the receiver is hugely complex. Instead of trying to mitigate frequency selective fading as a whole (which occurs when a huge bandwidth is allocated for the data transmission over a frequency selective fading channel), OFDM mitigates the problem by converting the entire frequency selective fading channel into small flat fading channels (as seen by the individual subcarriers). Flat fading is easier to combat (compared to frequency selective fading) by employing simple error correction and equalization schemes.

Difference between FDM and OFDM:
OFDM is a special case of FDM ( Frequency Division Multiplexing). In FDM, the given bandwidth is subdivided among a set of carriers. There is no relationship between the carrier frequencies in FDM. For example, consider that the given bandwidth has to be divided among 5 carriers (say a,b,c,d,e). There is no relationship between the subcarriers ; a,b,c,d and e can anything within the given bandwidth.

If the carriers are harmonics, say (b=2a,c=3a,d=4a,d=5a , integral multiple of fundamental component a ) then they become orthogonal. This is a special case of FDM, which is called OFDM (as implied by the word - 'orthogonal' in OFDM).

Eb/N0 Vs BER for BPSK over Rayleigh Channel and AWGN Channel

2011-05-06T08:46:00.002+01:00

The phenomenon of Rayleigh Flat fading and its simulation using Clarke's model and Young's model were discussed in the previous posts. The performance (Eb/N0 Vs BER) of BPSK modulation (with coherent detection) over Rayleigh Fading channel and its comparison over AWGN channel is discussed in this post.

We first investigate the non-coherent detection of BPSK over Rayleigh Fading channel and then we move on to the coherent detection. For both the cases, we consider a simple flat fading Rayleigh channel (modeled as a - single tap filter - with complex impulse response - h). The channel also adds AWGN noise to the signal samples after it suffers from Rayleigh Fading.

The received signal y can be represented as

where n is the noise contributed by AWGN which is Gaussian distributed with zero mean and unit variance and h is the Rayleigh Fading response with zero mean and unit variance. (For a simple AWGN channel without Rayleigh Fading the received signal is represented as y=x+n)

Non-Coherent Detection:

In non-coherent detection, prior knowledge of the channel impulse response ("h" in this case) is not know at the receiver. Consider the BPSK signaling scheme with x=+/- a being transmitted over such a channel as described above. This signaling scheme fails completely (in non coherent detection scheme), even in the absence of noise, since the phase of the received signal y is uniformly distributed between 0 and 2pi regardless of whether x[m]=+a or x[m]=-a is transmitted. So the non coherent detection of the BPSK signaling is not a suitable method of detection especially in a Fading environment.

Coherent Detection:
In coherent detection, the receiver has sufficient knowledge about the channel impulse response.Techniques like pilot transmissions are used to estimate the channel impulse response before the actual data transmission could begin. Lets consider that the impulse response estimate at receiver is perfect and accurate.
The transmitted symbols ('x') can be obtained from the received signal ('y') by the process of equalization as given below.

here z is still an AWGN noise except for the scaling factor 1/h. Now the detection of x can be performed in a
manner similar to the detection in AWGN channels.

The input binary bits to the BPSK modulation system are detected as

Theoretical BER:

The theoretical BER for BPSK modulation scheme over Rayleigh fading channel (with AWGN noise) is given by

The theoretical BER for BPSK modulation scheme over an AWGN channel is given by

Simulation : (click here to download Matlab code)

The following model is used for the simulation of BPSK over Rayleigh Fading channel and its comparison with AWGN channel

The Simulated and theoretical performance curves (Eb/N0 Vs BER) for BPSK modulation over Rayleigh Fading channel and the AWGN is given below.

Recommended Books:

Simulation of Rayleigh Fading ( Clarke's Model - sum of sinusoids method)

2011-05-02T07:22:00.001+01:00

A multipath fading channel can be modeled as a FIR (Finite Impulse Response) filter with the following impulse response.

where h(τ,t) is the time varying impulse response of the multipath fading channel having L multipaths
and h_i(t) and τ_i(t) denote the time varying complex gain and excess delay of the i-th path.

The above mentioned impulse response can be implemented as a FIR filter as shown below

The channel under consideration can be modeled as a multipath fading channel in which the impulse response may follow distributions like Rayleigh distribution ( in which there is no Line of Sight (LOS) ray between transmitter and receiver) or as Rician distribution ( dominant LOS path exist between transmitter and receiver), Nagami distribution, Weibull distribution etc.

Different methods of simulation techniques were proposed to simulate/model multipath channels. Some of the models include clarke's reference model, Jake's model, filteres gaussian noise model etc.

A Rayleigh fading channel (flat fading channel) is considered in this text.For simplicity we fix the excess delays τ_i(t) in the above equation and we generate h_i(t) that follows Rayleigh distribution. In this simulation Clarke's Rayleigh fading model is used. This model is also called mathematical reference model and is commonly considered as a computationally inefficient model compared to Jake's Rayleigh Fading simulator.

Theory of Rayleigh Fading:
Lets denote the complex impulse response h(t) of the flat fading channel as follows :

where hI(t) and hQ(t) are zero mean gaussian distributed. Therefore the fading envelope is Rayleigh distributed and is given by

The probability density function (Rayleigh distribution) of the above mentioned amplitude response is given by

We will use the Clarke's Rayleigh Fading model (given below) and check the statistical properties of the random process generated by the model against the statistical properties of Rayleigh distribution (given above).

Clarke's Rayleigh Fading model:
The random process of flat Rayleigh fading with M multipaths can be simulated with the sum-of-sinusoid method described as

Simulation : click here to download Matlab code
1) The rayleigh fading model is implemented as a function in matlab with following parameters:
M=number of multipaths in the fading channel, N = number of samples to generate, fd=maximum Doppler spread in Hz, Ts = sampling period

function [h]=rayleighFading(M,N,fd,Ts)

% function to generate rayleigh Fading samples based on Clarke's model

% M = number of multipaths in the channel

% N = number of samples to generate

% fd = maximum Doppler frequency

% Ts = sampling period

% Author : Mathuranathan Viswanathan for http://gaussianwaves.blogspot.com

% Creative commons License : CC-BY-NC-SA

a=0;
b=2*pi;

alpha=a+(b-a)*rand(1,M); %uniformly distributed from 0 to 2 pi
beta=a+(b-a)*rand(1,M); %uniformly distributed from 0 to 2 pi
theta=a+(b-a)*rand(1,M); %uniformly distributed from 0 to 2 pi

m=1:M;
for n=1:N;
     x=cos(((2.*m-1)*pi+theta)/(4*M));
     h_re(n)=1/sqrt(M)*sum(cos(2*pi*fd*x*n'*Ts+alpha));
     h_im(n)=1/sqrt(M)*sum(sin(2*pi*fd*x*n'*Ts+beta));
end
h=h_re+j*h_im;
end

2)The above mentioned function is used to generate Rayleigh Fading samples with the following values for the function arguments. M=15; N=10^5; fd=100 Hz;Ts=0.0001 second;

Investigation of Statistical Properties of samples generated using Clarke's model:
3) Mean and Variance of the real and imaginary parts of generated samples are
Mean of real part ~=0
Mean of imag part ~=0
Variance of real part = 0.4989 ~=0.5
Variance of imag part = 0.4989 ~=0.5

The results implies that the mean of the real and imaginary parts are same and are equal to zero.
The variance of the real and imaginary parts are approximately equal to 0.5.

4)Next, the pdf of the real part of the simulated samples are plotted and compared against the pdf of Gaussian distribution (with mean=0 and variance =0.5)

5) The pdf of the generated Rayleigh fading samples are plotted and compared against pdf of Rayleigh distribution ( with variance = 1)

6) From 4) and 5) we confirm that the samples generated by Clarke's model follows Rayleigh distribution (with variance = 1) and the real and imaginary part of the samples follow Gaussian distribution (with mean=0 and variance =0.5).

7) The Magnitude and Phase response of the generated Rayleigh Fading samples are plotted here.

More posts on Rayleigh Fading:
Rayleigh Fading simulation using Young's model
Introduction to Rayleigh Fading

References:
Marvin K Simon et al,"Digital communication over Fading Channels",Wiley-IEEE Press,2004
Theodore S. Rappaport,"Wireless Communications: Principles and Practice", Prentice Hall, 2002

Recommended Books:

Correlative Coding - Modified Duobinary Signaling

2011-04-26T07:30:00.002+01:00

Modified Duobinary Signaling is an extension of duobinary signaling. Modified Duobinary signaling has the advantage of zero PSD at low frequencies (especially at DC ) which is suitable for channels with poor DC response. It correlates two symbols that are 2T time instants apart, whereas in Duobinary signaling , symbols that are 1T apart are correlated.

The general condition to achieve zero ISI is given by

As discussed in the previous post, in correlative coding , the requirement of zero ISI condition is relaxed as a controlled amount of ISI is introduced in the transmitted signal and is counteracted in the receiver side

In the case of modified duobinary signaling, the above equation is modified as

which states that the ISI is limited to two alternate samples. Here a controlled or "deterministic" amount of ISI is introduced and hence its effect can be removed upon signal detection at the receiver.

Modified Duobinary Signaling:
The following figure shows the modified duobinary signaling scheme (click to enlarge).

Encoding Process:
1) a_n = binary input bit; a_n ∈ {0,1}.
2) b_n = NRZ polar output of Level converter in the precoder and is given by,

where a_k is the precoded output (before level converter).

3) y_n can be represented as

Note that the samples b_n are uncorrelated ( i.e either +d for "1" or -d for "0" input). On the other-hand,the samples y_n are correlated ( i.e. there are three possible values +2d,0,-2d depending on ak and a_k-2). Meaning that the modified duobinary encoding correlates present sample a_k and the previous input sample a_k-2.

4) From the diagram,impulse response of the modified duobinary encoder is computed as

Decoding Process:
5) The receiver consists of a modified duobinary decoder and a postcoder (inverse of precoder).The modified duobinary decoder implements the following equation (which can be deduced from the equation given under step 3 (see above))

This equation indicates that the decoding process is prone to error propagation as the estimate of present sample relies on the estimate of previous sample. This error propagation is avoided by using a precoder before modified-duobinary encoder at the transmitter and a postcoder after the modified-duobinary decoder. The precoder ties the present sample and the sample that precedes the previous sample ( correlates these two samples) and the postcoder does the reverse process.

6) The entire process of modified-duobinary decoding and the postcoding can be combined together as one algorithm. The following decision rule is used for detecting the original modified-duobinary signal samples {a_n} from {y_n}

Matlab Code: - click here to download

Output:

Input Symbols to the system : (an)
1 0 1 1 0 0 1 1 1 1

Decoded Symbols : (estimated an)
1 0 1 1 0 0 1 1 1 1

Recommended Books:

Correlative coding - Duobinary Signaling

2011-04-21T04:01:00.002+01:00

The condition for zero ISI (Inter Symbol Interference) is

which states that when sampling a particular symbol (at time instant nT=0), the effect of all other symbols on the current sampled symbol is zero.

As discussed in the previous post, one of the practical ways to mitigate ISI is to use partial response signaling technique ( otherwise called as "correlative coding"). In partial response signaling, the requirement of zero ISI condition is relaxed as a controlled amount of ISI is introduced in the transmitted signal and is counteracted in the receiver side.

By relaxing the zero ISI condition, the above equation can be modified as,

which states that the ISI is limited to two adjacent samples. Here we introduce a controlled or "deterministic" amount of ISI and hence its effect can be removed upon signal detection at the receiver.

Duobinary Signaling:
The following figure shows the duobinary signaling scheme (click to enlarge).

Encoding Process:
1) a_n = binary input bit; a_n ∈ {0,1}.
2) b_n = NRZ polar output of Level converter in the precoder and is given by,

3) y_n can be represented as

Note that the samples b_n are uncorrelated ( i.e either +d for "1" or -d for "0" input). On the other-hand,the samples y_n are correlated ( i.e. there are three possible values +2d,0,-2d depending on a_n and a_n-1). Meaning that the duobinary encoding correlates present sample a_n and the previous input sample a_n-1.

4) From the diagram,impulse response of the duobinary encoder is computed as

Decoding Process:
5) The receiver consists of a duobinary decoder and a postcoder (inverse of precoder).The duobinary decoder implements the following equation (which can be deduced from the equation given under step 3 (see above))

This equation indicates that the decoding process is prone to error propagation as the estimate of present sample relies on the estimate of previous sample. This error propagation is avoided by using a precoder before duobinary encoder at the transmitter and a postcoder after the duobinary decoder. The precoder ties the present sample and previous sample ( correlates these two samples) and the postcoder does the reverse process.

6) The entire process of duobinary decoding and the postcoding can be combined together as one algorithm. The following decision rule is used for detecting the original duobinary signal samples {a_n} from {y_n}

Matlab Code: - click here to download

Output:

Input Symbols to the system : (an)
0 0 1 0 1 1 0 1 1 1

Decoded Symbols : (estimated an)
0 0 1 0 1 1 0 1 1 1

Recommended Books:

Introduction to controlled ISI (Inter Symbol Interference)

2011-04-19T11:51:00.005+01:00

Most of the communication channels (data storage devices,optical,wireless channels etc..) can be considered as band-limited linear filters. The channels can be modeled as having the following frequency response (over a given bandwidth W).

Here A(f) is amplitude response and θ(f) is the phase response of the channel over the given bandwidth W.
The envelope or group delay is defined as,

A channel is considered non-distorting (within the given bandwidth W occupied by the transmitted signal), when the amplitude response is constant and the phase response is a linear function frequency (within the given bandwidth W) (or the group delay is a constant).

Amplitude distortion occurs when the amplitude response is no longer constant.Delay/phase distortion occurs when the phase response is not a linear function of frequency or the envelope/group delay is not a constant.

A non-ideal channel frequency response is caused by amplitude and phase distortion. When a succession of pulses transferred through a non-ideal channel, at a rate 2W symbols/second (R=1/T=2W -'Nyquist Rate'), gets distorted to a point that they are no longer distinguishable from each other. This is called Inter Symbol Interference (ISI), meaning that a symbol transmitted across a non-ideal channel will be affected by the other symbols. Pulse shaping filters are generally employed at the transmitter and receiver to match the spectral characteristics of the signal with that of the channel.

Generally following strategies may be employed to mitigate ISI

1) Use ideal rectangular pulse shaping filters to achieve zero ISI:
Maximum transferable data rate that is possible with zero ISI is R=1/T=2W symbols/second, provided ideal rectangular transmit and receiver pulse shaping filters are used. Ideal rectangular transmit and receive filters are practically unrealizable. So this option is not viable to achieve zero ISI.

See also: Implementation of rectangular pulse shaping filters

2) Relax the condition of transmitting at maximum rate R=1/T=2W, to achieve zero ISI:
If the transmission rate is reduced below 2W (i.e R=1/T<2W), then it is possible to implement practically realizable filters. Raised Cosine and Square Root Raised cosine filters are generally used to achieve zero ISI, if the transmission data rate is reduced below 2W. The signals generated using this method are called full response signals.

See also: Implementation of Square Root Raised Cosine pulse shaping filters
Implementation of Raised Cosine pulse shaping filters

3) Relax the condition of zero ISI and transfer at Nyquist Rate (R=1/T=2W):
In this method, we relax the condition of achieving zero ISI , so that the data can be transferred at maximum possible rate ( R=1/T=2W).Instead of achieving zero ISI, this method introduces controlled amount of ISI in the transmitted signal and counteracts it upon receiving it. The transmit filter is designed to introduce 'deterministic' or 'controlled' amount of ISI and is counteracted in the receiver side. Methods like duobinary signaling, modified duobinary signaling are employed under this category. The resulting signals are called partial response signals which are transmitted at Nyquist rate of 2W symbols/second.

See also: Implementation of duobinary signaling
Implementation of modified duobinary signaling

Recommended Books

Square Root Raised Cosine Filter (Matched/split filter implementation)

2011-04-07T11:02:00.001+01:00

An implementation of a rectangular pulse shaping filter in transmitter and receiver side was discussed in the previous post.

The code was modified a bit to implement a Square Root Raised Cosine Filter (SRRC). Instead of using a single Raised Cosine filter at the transmitter, a square root raised cosine filter is used at both transmitter and receiver. This is to achieve matched filtering (maximizes SNR of the system). SRRC and RRC filters are used to mitigate inter symbol interference (ISI).

The frequency response of the SRRC filter is the square root of the frequency response of the raised-cosine filter. The impulse response of the SRRC filter (according to 3GPP standard[1]) is given by

with following values^[2] at t=0 and t/T=±1/4α .

There are other forms of equations to implement SRRC.

Matlab function to implement Square Root Raised Cosine Filter:

function [response]=srrc(os_factor,roll_off)

%os_factor=oversampling factor

a=roll_off;
t=-4:1/os_factor:4; %Limiting the response to -4T to 4T
%This can be increased or decreased according to the requirement
p=zeros(1,length(t));
    for i=1:1:length(t)
        if t(i)==0
            p(i)= (1-a)+4*a/pi;
        else if t(i)==1/(4*a) || t(i)==-1/(4*a)
           p(i)=a/sqrt(2)*((1+2/pi)*sin(pi/(4*a)+(1-2/pi)*cos(pi/(4*a))));
              else
                p(i) = (sin(pi*t(i)*(1-a))+4*a*t(i).*cos(pi*t(i)*(1+a)))./(pi*t(i).*(1-(4*a*t(i)).^2));
             end
        end
    end
    response=p./sqrt(sum(p.^2)); %Normalization to unit energy
end

Simulation of Square Root Raised Cosine Filter in Matlab:

Define oversampling factor (to simulated digital to analog conversion in the transmitter), input bits to be transmitted and the impulse response of the desired SRRC filter (as shown in the equation/function above).
The SRRC filter's impulse response is computed for t=-4T to 4T duration. This duration can be changed to any appropriate desired value.

clc;clear;close all;
overSampling_Factor=8;
Input_bit = [1]; %Bits to be transmitted
Input_bit_os=upsample(Input_bit,overSampling_Factor); %oversampling
alpha=0.1; % roll-off factor of Root Raised Cosine Filter
pt = srrc(overSampling_Factor,alpha); % impulse response of SRRC filter

Convolve the impulse response of the root raised cosine filter with the oversampled version of input bits and plot the resulting response. This is the response of the filter for the input bits at the transmitter side.

output_of_srrc_filter = conv(Input_bit_os,pt);
stem(output_of_srrc_filter);
title('Response of SRRC Filter at Tx side')
xlabel('Samples')
ylabel('Amplitude')

At the receiver side the received signal ("output_of_srrc_filter") is convolved with a matched filter (which is same as the transmitting filter).Note that no noise is added to the incoming signal. If you wish you could add some noise to the signal before performing the matched filter operation.

%Receiver side; Using a matched filter (that is matched to the SRRC pulse

%in the transmitter)

y = conv(output_of_srrc_filter,pt);
stem(y);
title('Matched filter (SRRC) response at Rx side')
xlabel('Samples')
ylabel('Amplitude')

Note that the peak value of the matched filter response occurs at 65th sample (this number will vary depending on the oversampling rate). Noting down the sample number of the peak value for a single bit input at the transmitter is important in downsampling operation (equivalent to recovering the input bits by sampling using a Digital to Analog Converter). A small snipped of code is given below, to extract the recovered bits out of the matched filter response.

midSample=length(-4:1/overSampling_Factor:4);
y_truncated=y(midSample-1:end); %Remove unwanted portions(first few samples till the peak value)

%Now the first sample contains the peak value of the response. From here

%the samples are extracted depending on the oversampling factor

y_down = downsample(y_truncated,overSampling_Factor,1); %here offset=1 means starting from 1st sample %retain every 8th sample

stem(y_down);
title('Downsampled output (ADC conversion and Sampling)')
xlabel('Samples')
ylabel('Amplitude')

Following simulation results are shown for a series of input bits: Input_bit = [1 -1 1 -1 1 1 -1 -1 -1 1 -1 1 1 1]. It is interesting to note how the effect of inter symbol interference being mitigated by the Square Root Raised Cosine Filters implemented at both the transmitter and receiver. The roll-factor of the SRRC filter can be adjusted to see the effects of ISI on the recovered samples.

(Click the picture to enlarge)

References:

[1] 3GPP TS 25.104 V6.8.0 (2004-12) (URL for all 3GPP Standards - click here)
[2]Michael Joost,"Theory of Root-Raised Cosine Filter",Krefeld DE,December 2010

Recommended Books

Oversampling, ADC - DAC Conversion,pulse shaping and Matched Filter

2011-04-07T06:01:00.000+01:00

In order to simulate a communication system in Matlab one needs to understand the concept of Oversampling (upsampling) and undersampling (downsampling).

Digital communication systems involves conversion of digital data to analog form (with some modulation,coding stuffs etc..) at the transmitter side. The data transmitted through a channel is in analog form and it needs to be converted back to digital domain by the process of Digital to Analog Conversion (DAC). Once the data is in digital domain, the receiver is at ease to process it to recover the transmitted data.

In Matlab one needs to emulate this situation where in the input data has to be represented in analog form in the transmitter and back to digital form at the receiver. Matlab is just a piece of software that runs on a digital machine (computer terminals), so it cannot really produce data in analog domain. To circumvent this problem, one has to use a process called oversampling.

Lets consider a digital communication system (shown below), where in a train of input data (In = +/-1) are shaped by a pulse shaping filter and modulated by a carrier Fc. This is a simple system implementing BPSK modulation. We will use rectangular pulse shaping filter (p(t)) for implementation in Matlab.

For illustration we try to transmit only one bit (In=1) through the system. To represent digital to analog conversion ("n" domain to "t" domain), we use oversampling by 4. Oversampling is a process in which additional "oversampling factor-1" zeros are inserted in to the digital data.

For example, oversampling In=1 by a factor of 4 gives [1 0 0 0]. Similarly, oversampling the input [1 0 1 1] by a factor of 3 gives [1 0 0 0 0 0 1 0 0 1 0 0].

overSampling_Factor=4;
Input_bit = 1;
Input_bit_os=upsample(Input_bit,overSampling_Factor); %oversampling by 4

%This is equivalent to performing digital to analog conversion in Tx side

The oversampled version of input is then passed through a pulse shaping filter (it can be a rectangular filter, root cosine filter, etc..,). A rectangular filter is implemented in the following code.

%Normalizing the pulse shape to have unit energy

pt = [ones(1,overSampling_Factor) 0 0 0 0 0 0]/sqrt(overSampling_Factor);

% impulse response of a rectangular pulse

%convolving the oversampled input with rectangular pulse

%The output of the convolution operation will be in the transmitter side

output_of_rect_filter = conv(Input_bit_os,pt);

stem(output_of_rect_filter);
title('Output of Rectangular Filter at Tx side')
xlabel('Samples')
ylabel('Amplitude')

At the receiver side, we might want to recover the transmitter samples. To do this we need to used a matched filter at the receiver side (to maximize the SNR of the system). The matched filter is exactly matched to the pulse shaping filter used in the transmitter side.

%Receiver side; Using a matched filter (that is matched to the rect pulse

%in the transmitter)

y = conv(output_of_rect_filter,pt);
stem(y)
title('Matched filter output at Rx side')
xlabel('Samples')
ylabel('Amplitude')

Once the received signal is passed through the matched filter it has to be converted back to digital domain ("t" to "n" domain). This is achieved by implementing a downsampling process. Note that the value "1" is at 4 th sample, this is where we have to sample in the receiver side.

%Downsampling by 4, since the actual value of the output is shifted to 4 th
%sample
y_down = downsample(y,overSampling_Factor,overSampling_Factor-1);

stem(y_down)
title('Downsampled output (ADC conversion and Sampling)')
xlabel('Samples')
%ylabel('Amplitude')

The above chart illustrates the scenario in which the input bit is received properly via a matched filter arrangement.Note that there are some trailing zeros in the downsampled signal (which is the effect of padding additional zeros in the pulse shaping function). The trailing zeros can be removed, if the padded zeros are stripped off from the pulse shaping function ( Matlab Code: pt = [ones(1,overSampling_Factor) ]/sqrt(overSampling_Factor);)

Output of the received signal (after matched filter) for an input [1 0 1 1] is shown below.

Recommended Products/books on Matlab:

BER Vs Eb/N0 for 8-PSK modulation over AWGN

2011-03-25T09:26:00.002Z

BER Vs Eb/No for QPSK and BER were discussed in the previous posts.Click below to see the posts if you have missed them earlier.
(BER Vs Eb/N0 for QPSK modulation over AWGN),
(BER Vs Eb/N0 for BPSK modulation over AWGN)

BER Vs Eb/No for 8-PSK is discussed here. A baseband simulation of 8-PSK is demonstrated below.

Simulation Model: (click to enlarge the figure)

Transmitter:

A complex baseband M-PSK signal is represented by

where Es is the symbol energy of each M-PSK Symbol, T is the bit period, θn is the phase shift for each of the symbol. For 8-PSK, M=8.

In 8PSK modulator, every 3 input data bits are mapped to one complex-valued 8PSK symbol (where as in QPSK 2 input bits are mapped to one complex-valued QPSK symbol). Each symbol is gray coded prior to mapping them to the 8-PSK constellation (this is to allow only one bit change in the constellation, there by avoiding abrupt phase changes in the constellation which may lead to addition errors in the receiver) . In 8PSK there are eight possible symbols that can be transmitted. We will use eight different phases ( 22.5°, 67.5°, 112.5°, 157.5°, 202.5°, 247.5°, 292.5°, 337.5° ) to represent each 8-PSK symbol. Each symbol is 45° away from adjacent symbols (see the constellation diagram below).

The 8PSK symbols are transmitted over an AWGN channel that adds Gaussian Noise to the transmitted symbol.

Receiver:

The demodulation process is performed by finding one of the eight constellation symbols that has the smallest Euclidean distance with received symbol. Some may prefer to use a sector based detection approach in which the received signal constellation is divided into sectors.For 8PSK, the entire constellation can be equally divided into 8 sectors, with each angle sector corresponding to one of the 8 constellation symbols. The demodulation process can be performed by using the angle sectors. For example, the angle sector defined by [0,π/4] corresponds to the 8PSK symbol S₁. Instead of this approach, this simulation is based on minimum Euclidean Distance. The demodulator computes the Euclidean Distance of each received symbol against all possible ideal symbols. Then it pin points an ideal symbol which gave the minimum Euclidean Distance.

Euclidean Distance is given by

Where y=received signal (with noise), s=points on ideal constellation and the subscripts R and I are for their real and imaginary part respectively.
The demodulator calculates the Euclidean distance of a received symbol against all possible symbols in the constellation and picks a symbol from the ideal constellation that gives the minimum Euclidean Distance.

Finally, the demodulated symbols are converted back to binary symbols and the BER (Bit Error Rate) is calculated.

Matlab Program: (Click here to download the entire package)

%   8PSK Simulation
%   Author: Mathuranathan Viswanathan (http://gaussianwaves.blogspot.com)
%   License : creative commons : Attribution-NonCommercial-ShareAlike 3.0
%   Unported
clc;
clear;
close all;

N=100000; %Number of data bits to send over the channel
EbN0dB=-6:2:12;

N=N+rem((3-rem(N,3)),3); %add additional bits to the data to make the length multiple of 3
% one 8-PSK symbol contains 3 binary bits
x=rand(1,N)>=0.5;%Generate random 1's and 0's as data;

%Club 3 bits together and gray code it individually
inputSymBin=reshape(x,3,N/3)';
g=bin2gray(inputSymBin);

%Convert each Gray coded symbols to decimal value this is to ease the process
%of mapping based on array index
b=bin2dec(num2str(g,'%-1d'))';

%8-PSK Constellation Mapper
%8-PSK mapping Table
map=[1.414 0.707;0.707 1.414;-0.707 1.414;-1.414 0.707;-1.414 -0.707;-0.707 -1.414;0.707 -1.414;1.414 -0.707];
s=map(b(:)+1,1)+1i*map(b(:)+1,2);

%Simulation for each Eb/N0 value
M=8; %Number of Constellation points M=2^k for 8-PSK k=3
Rm=log2(M);   %Rm=log2(M) for 8-PSK M=8
Rc=1;   %Rc = code rate for a coded system. Since no coding is used Rc=1

simulatedBER = zeros(1,length(EbN0dB));
theoreticalBER = zeros(1,length(EbN0dB));
count=1;
figure;
for i=EbN0dB
    %-------------------------------------------
    %Channel Noise for various Eb/N0
    %-------------------------------------------
    %Adding noise with variance according to the required Eb/N0
    EbN0 = 10.^(i/10); %Converting Eb/N0 dB value to linear scale
    noiseSigma = sqrt(2)*sqrt(1./(2*Rm*Rc*EbN0)); %Standard deviation for AWGN Noise
    %Creating a complex noise for adding with 8-PSK s signal
    %Noise is complex since 8-PSK is in complex representation
    n = noiseSigma*(randn(1,length(s))+1i*randn(1,length(s)))';
    y = s + n;

    plot(real(y),imag(y),'r*');hold on;
    plot(real(s),imag(s),'ko','MarkerFaceColor','g','MarkerSize',7);hold off;
    title(['Constellation plots - ideal 8-PAM (green) Vs Noisy y signal for EbN0dB =',num2str(i),' dB']);
    pause;

    %Demodulation
    %Find the signal points from MAP table using minimum Euclidean distance
    demodSymbols = zeros(1,length(y));
    for j=1:length(y)
        [minVal,minindex]=min(sqrt((real(y(j))-map(:,1)).^2+(imag(y(j))-map(:,2)).^2));
        demodSymbols(j)=minindex-1;
    end

    demodBits=dec2bin(demodSymbols)-'0'; %Dec to binary vector

    xBar=gray2bin(demodBits)'; %gray to binary
    xBar=xBar(:)';

    bitErrors=sum(sum(xor(x,xBar)));
    simulatedBER(count) = log10(bitErrors/N);
    theoreticalBER(count) = log10(1/3*erfc(sqrt(EbN0*3)*sin(pi/8)));

    count=count+1;
end
figure;
plot(EbN0dB,theoreticalBER,'r-*');hold on;
plot(EbN0dB,simulatedBER,'k-o');

title('BER Vs Eb/N0 (dB) for 8-PSK');legend('Theoretical','Simulated');grid on;
xlabel('Eb/N0 dB');
ylabel('BER - Bit Error Rate');
grid on;

Output:

Recommended Books:

Walsh Hadamard Code - Matlab Simulation

2011-03-03T02:47:00.006Z

The following is a function to generate a Walsh Hadamard Matrix of given codeword size. The codeword size has to be a power of 2.

function [H]=generateHadamardMatrix(codeSize)

%[H]=generateHadamardMatrix(codeSize);

% Function to generate Walsh-Hadamard Matrix where "codeSize" is the code

% length of walsh code. The first matrix gives us two codes; 00, 01. The second

% matrix gives: 0000, 0101, 0011, 0110 and so on

% Author: Mathuranathan Viswanathan for http://gaussianwaves.blogspot.com

% License: Creative Commons: Attribution-NonCommercial-ShareAlike 3.0

% Unported

%codeSize=64; %For testing only

N=2;
H=[0 0 ; 0 1];
if bitand(codeSize,codeSize-1)==0
while(N~=codeSize)
       N=N*2;
       H=repmat(H,[2,2]);
       [m,n]=size(H);

%Invert the matrix located at the bottom right hand corner

      for i=m/2+1:m,
      for j=n/2+1:n,
                H(i,j)=~H(i,j);
         end
     end
end
else
disp('INVALID CODE SIZE:The code size must be a power of 2');
end

Example:

To Generate Walsh Codes used in IS-95 (which utilizes 64 Walsh codes of size 64 bits each, use : [H]=generateHadamardMatrix(64)
This will generate 64 Walsh Codes of length 64-bits (for each code).

Test Program: Click here to download the entire package
Also given below is a program to test the cross-correlation and auto-correlation of Walsh code. A set of 8-Walsh codes are used for this purpose.

% Matlab Program to test Walsh Hadamard Codes and to test its orthogonality

% Plots cross-correlation and auto correlation of Walsh Hadamard Codes

% Author: Mathuranathan Viswanathan for http://gaussianwaves.blogspot.com

% License: Creative Commons: Attribution-NonCommercial-ShareAlike 3.0

% Unported

codeSize=8;
[H]=generateHadamardMatrix(codeSize);

%-----------------------------------------------------------

%Cross-Correlation of Walsh Code 1 with rest of Walsh Codes

h = zeros(1, codeSize-1); %For dynamic Legends
s = cell(1, codeSize-1); %For dynamic Legends
for rows=2:codeSize
[crossCorrelation,lags]=crossCorr(H(1,:),H(rows,:));
h(rows-1)=plot(lags,crossCorrelation);
s{rows-1} = sprintf('Walsh Code Sequence #-%d', rows);
hold all;
end

%Dynamic Legends

% Select the plots to include in the legend

index = 1:codeSize-1;

% Create legend for the selected plots

legend(h(index),s{index});
title('Cross Correlation of Walsh Code 1 with the rest of the Walsh Codes');
ylabel('Cross Correlation');
xlabel('Lags');

%-----------------------------------------------------------

%AutoCorrelation of Walsh Code - 1

autoCorr2(H(2,:),8,2,1);

Click here to download the entire package

Output:

From the output plots below, it can be ascertained that the Walsh codes has excellent cross-correlation property and poor autocorrelation property. Excellent cross-correlation property (zero cross-correlation) implies orthogonality, which makes it suitable for CDMA applications

More on Walsh Codes and CDMA:
Codes used in CDMA (Long codes, Short Codes and Walsh Codes)
Generation of Gold Codes and their Cross-correlation
Preferred Pairs m-sequence generation for Gold Codes
Maximum Length Sequences
Spread Spectrum Communications - Intro

Recommended Books on CDMA and spread spectrum techniques

Codes used in CDMA

2011-02-28T05:36:00.002Z

IS-95 CDMA Standard uses three types of codes namely 1) Long code 2) Short code and 3) Walsh Hadamard codes.

Long Code:

Long codes run at 1.2288 Mb/s and are 42 bits longs (created using a 42 bits LFSR register). It takes approx 41.2 days to repeat a long code at this rate. It is used for both encryption and spreading. Encryption is achieved by using a mask called Long Code mask which is a created using a 64-bit authentication key called A-key (assigned by CAVE protocol) and Electronic Serial Number ( ESN - assigned each user based on the mobile number). The Long code changes each time a new connection is created.

Short Code:

Short code is a m-sequence of lenght 2¹⁵-1 (created using a 15 bit LFSR register) and is used for synchronization purpose in the forward as well as reverse links. The short code is also used to identify cell/base station connection in the forward link. It repeats approx 75 times in 2 seconds. Each base station is assigned a cyclically shifted version of same short code sequence to differentiate the base stations.This is also called PN offset in CDMA jargon. Since the cyclically shifted versions of a same m-sequence offer poor correlation, it is easier to differentiate between different base station links.

During the initial call setup stage, a mobile phone tries to find a base station (in 2 seconds max allowed time), if it find a base station, the mobile phone is validated using a database by the base station and is assigned a PN Short code sequence. This PN short code sequence uniquely identifies the connection between the particular base station and the mobile devices served under that base station.
In reality two short code sequences are used one for I channel and another for Q channel (used in spreading and de-spreading).

Walsh Hadamard Code:

CDMA used another type of code called Walsh Hadamard Code. In IS-95 CDMA, 64 Walsh codes are used per base station. This enables to create 64 separate channels per base stations (i.e. a base station can handle maximum 64 unique users at a given time).In CDMA-2000 standard, 256 Walsh codes are used to handle maximum 256 unique users under a base.

Walsh codes are created using Hadamard Matrix and transform. The codes under a family of Walsh codes, posses a beautiful property of being orthogonal to each other (what else do we want to identify/accommodate users ?).

The first matrix in a Hadamard transform is

The next matrices are formed iteratively using

For example H2 will be

Each row of a Hadamard matrix represent a unique Walsh code and all the Walsh codes in a given matrices are orthogonal. The length of the row of the matrix ( number of columns otherwise) is the code-length of the Walsh codes. To get a 64-Walsh code matrix we need to transform the matrices till H8 (this matrix contains 64 rows - representing 64 walsh codes and each code is of 64 bits length).

Walsh codes posses excellent cross-correlation property ( cross correlation of one Walsh code with another is always zero) therefore possess excellent orthogonality property. The auto-correlation property of Walsh code is very poor and so it is used only in synchronous CDMA networks, which maintains a synchronizing mechanism to identify the starting of the codeword.

Actually in IS-95, out of the 64 available Walsh codes, Walsh code 0 is reserved for pilot channel, 1 to 7 are assigned for synch channel and paging channels and the remaining 8-63 are assigned for users (traffic channel).

Window Functions - An Analysis

2011-02-10T08:57:00.001Z

Effect of limited measurement time and it consequence - Spectral Leakage were discussed in the previous post. Lets discuss the application of window techniques that has the potential of reducing the spectral leakage.

The usage of windows (like Hann,Hamming,etc.., except for rectangular window) smoothens the observed signal over the edges of measured time interval and gradually reduces it to zero. This prevents glitches in the "assumed" signal repetitively reconstructed by the Fourier Transform. This phenomena reduces the spectral leakage.

The process of multiplying the measured signal with a smoothly ending function (that gradually reduces to zero) is called "windowing" technique and the function is called "window" function. Frequency domain characteristics of window functions are often analyzed to determine their suitability for a specific application/system limitations. Spectral characteristics of window functions are used to underpin the answer to the most commonly occurring dilemma of making a choice between "detection" and "resolution".

"Detection" means detecting a desired signal in the presence of broadband noise. Resolution refers to the ability of distinguishing narrowband spectral components. Available are plethora of parameters that would help us to select the window functions which would suit our needs.

Commonly used Parameters for Window Selection:^[1][2][3]

Main Lobe Width:
Main Lobe Width is directly linked to "Frequency resolution". Narrower the main lobe, more will be its frequency resolution (Ability of selecting finest frequency components). Main Lobe Width is the measure of the main lobe width at -3dB or -6dB below main lobe peak. When the main lobe width decreases, the remaining energy spreads out to side lobes thereby increasing spectral leakage/decreasing amplitude accuracy (decreasing "detection" ability). A compromise is needed to strike a balance between detection and resolution that would suit the application at hand.

Side Lobe Level:
Side lobes occur at either side of the main lobe and they approach zero at integral multiples of Fs/N (where Fs=Sampling Frequency; N is the length of N-point FFT/DFT). Side lobes directly affect the extent to which adjacent frequency components leak into adjacent frequency bins. Side Lobe Level is usually measure for that side lobe which has the maximum peak (compared to other side lobes). It is measured in decibels relative to the peak of the main lobe. Lesser the side lobe level of a window, more will be its detection ability.The side lobe roll-off rate is the asymptotic decay rate in decibels per decade of frequency of the peaks of the side lobes.(Note: Main Lobe Width and Side Lobe levels can be measured by visual inspection of the frequency response of the windows)

For an "Ideal" Window :
1) Main lobe width is very small (High Frequency Resolution)
2) Side Lobe level is very high (Good noise suppression - High detection ability)
3) Side Lobes fall-offs rapidly

Normalized Equivalent Noise Bandwidth:^[2]
The windowing of a signal in time also affects the resolution bandwidth that can be achieved. When calculating the spectrum of a signal segment the resolution bandwidth achieved with a window is always lower (worse) than the resolution bandwidth achieved without a window. The normalized equivalent noise bandwidth (NENBW) of a window is one measure of how much it reduces the resolution bandwidth that can be achieved. To compensate for this a longer signal segment should be processed.Equivalent noise bandwidth (ENBW) compares the window to an ideal, rectangular filter. It is the equivalent width of a rectangular filter that passes the same amount of white noise as the window. The normalized ENBW (NENBW) is the ENBW multiplied by the time duration of the of the signal being windowed. NENBW values for windows available in SpectrumAnalyzerResBW are listed in the following table.For our computation we assume that the time duration is of unit duration ,therefore ENBW and NENBW will be equal.

Calculation of Equivalent Noise Bandwidth:
Equivalent Noise Bandwidth is calculated as a ratio of inherent power gain over coherent power gain as given below. This equation is particularly useful for simulations as we do not have to bother about the transfer functions of each window functions in frequency domain.

List of Window Functions:^[3]

Matlab Code: (click here to download)

A Matlab Code for simulating the Window functions and its response is given below. The Main lobe level, Side lobe level ENBW and other parameters are measured based on the plots obtained and are tabulated below.

% Simulation of Window Functions and its characteristics
%Author: Mathuranathan Viswanathan (http://gaussianwaves.blogspot.com)
% License : creative commons : Attribution-NonCommercial-ShareAlike 3.0
% Unported
clear;
clc;
N=64;
n=0:N-1;
%Window Functions
w_rec=ones(1,N);
w_hann=0.5*(1-cos(2*pi*n/(N-1)));
w_hamming=0.54-0.46*cos(2*pi*n/(N-1));
w_bartlett=2/N*(N/2-abs(n-(N-1)/2));

alpha=0.16;a0=(1-alpha)/2;a1=1/2;a2=alpha/2;
w_blackman=a0-a1*cos(2*pi*n/(N-1))+a2*cos(4*pi*n/(N-1));

alpha=3;
num=besseli(0,pi*alpha*sqrt(1-(2*n/(N-1)-1).^2));
den=besseli(0,pi*alpha);
w_kaiser=num/den;

alpha=0.5;
w_gaussian=exp(-1/2*((n-(N-1)/2)./(alpha*(N-1)/2)).^2);

a0=1;a1=1.93;a2=1.29;a3=0.388;a4=0.032;
w_flattop=a0-a1*cos(2*pi*n/(N-1))+a2*cos(4*pi*n/(N-1))-a3*cos(6*pi*n/(N-1))+a4*cos(8*pi*n/(N-1));
w_flattop=w_flattop/max(w_flattop);

%Frequency Response of Window Functions
[W_REC,f]=freqz(w_rec,1);
W_REC=W_REC/(max(abs(W_REC)));
[W_HANN,f]=freqz(w_hann,1);
W_HANN=W_HANN/(max(abs(W_HANN)));
[W_HAMMING,f]=freqz(w_hamming,1);
W_HAMMING=W_HAMMING/(max(abs(W_HAMMING)));
[W_BARTLETT,f]=freqz(w_bartlett,1);
W_BARTLETT=W_BARTLETT/(max(abs(W_BARTLETT)));
[W_BLACKMAN,f]=freqz(w_blackman,1);
W_BLACKMAN=W_BLACKMAN/(max(abs(W_BLACKMAN)));
[W_KAISER,f]=freqz(w_kaiser,1);
W_KAISER=W_KAISER/(max(abs(W_KAISER)));
[W_GAUSSIAN,f]=freqz(w_gaussian,1);
W_GAUSSIAN=W_GAUSSIAN/(max(abs(W_GAUSSIAN)));
[W_FLATTOP,f]=freqz(w_flattop,1);
W_FLATTOP=W_FLATTOP/(max(abs(W_FLATTOP)));

%Plotting
figure;
handle=plot(n,[w_rec;w_hann;w_hamming;w_bartlett;w_blackman;w_kaiser;w_gaussian]);
set(handle,'linewidth',1.5);hold on;
handle=plot(n,[w_flattop],'b-*');
set(handle,'linewidth',1.5);axis([0 N-1 -0.2 1.2]);grid on;
legend('Rectangular','Hanning','Hamming','Bartlett','Blackman','Kaiser','Gaussian','Flat-top',0);
title('Window Functions in Time domain - w(n)');
%xlabel('Samples');ylabel('Amplitude');

figure;
handle=plot(20*log10(abs(([W_REC W_HANN W_HAMMING W_BARTLETT W_BLACKMAN W_KAISER W_GAUSSIAN]))));
set(handle,'linewidth',1.5);hold on;
handle=plot(20*log10(abs(([W_FLATTOP]))),'b-*');
set(handle,'linewidth',1.5);axis([0 length(W_REC)-1 -200 1]);grid on;
legend('Rectangular','Hanning','Hamming','Bartlett','Blackman','Kaiser','Gaussian','Flat-top',0);
title('Frequency Response of Window function- W(K)');
xlabel('Frequency Bins (K)');%ylabel('Magnitude');

%Equivalent Noise Bandwidth Calculations
ENBW_Rectangular=(1/N*sum(w_rec.^2))/((1/N*sum(w_rec)).^2)
ENBW_Hanning=(1/N*sum(w_hann.^2))/((1/N*sum(w_hann)).^2)
ENBW_Hamming=(1/N*sum(w_hamming.^2))/((1/N*sum(w_hamming)).^2)
ENBW_Bartlett=(1/N*sum(w_bartlett.^2))/((1/N*sum(w_bartlett)).^2)
ENBW_Blackman=(1/N*sum(w_blackman.^2))/((1/N*sum(w_blackman)).^2)
ENBW_Kaiser=(1/N*sum(w_kaiser.^2))/((1/N*sum(w_kaiser)).^2)
ENBW_Gaussian=(1/N*sum(w_gaussian.^2))/((1/N*sum(w_gaussian)).^2)
ENBW_Flattop=(1/N*sum(w_flattop.^2))/((1/N*sum(w_flattop)).^2)

%Other functions to try
%Bartlett-Hann Window
%a0=0.62;a1=0.48;a2=0.38;
%w_bartlett_hann=a0-a1*abs(n/(N-1) - 1/2)-a2*cos(2*pi*n/(N-1));
%W_BARTLETT_HANN=freqz(w_bartlett_hann,1,f);
%W_BARTLETT_HANN=W_BARTLETT_HANN/(max(abs(W_BARTLETT_HANN)));

Simulation Results and Tabulation of "Measured" Window Parameters:

References:

[1] National Instruments, "Characteristics of Different Smoothing Windows" ,LabView 8.2 Help,Aug 2006

[2] Agilent Technologies,"Spectrum analyzer with resolution bandwidth setting",Advanced Design System Documentation 2008.
[3] Wikipedia, "Window Function" ,Wikimedia Foundation Inc.

Additional Resources:

[1] If you are particularly interested in calculating the Equivalent Noise Bandwidth of Analog Filters refer - Thomas J. Karras,"Equivalent Noise Bandwidth Analysis from Transfer Functions",NASA Technical Note,NASA TN D-2842

[2] Strategies for Choosing Windows : Michael Cerna and Audrey F. Harvey,"The Fundamentals of FFT-Based Signal Analysis and Measurement",National Instruments Application Notes 041.

FFT and Spectral Leakage

2011-01-20T10:34:00.001Z

Introduction:

Frequency Transform is used to study a signal's frequency domain characteristics. When using FFT to study the frequency domain characteristics of a signal, there are two limits : 1) The detectability of a small signal in the presence of a larger one ; 2) frequency resolution - which distinguishes two different frequencies.

In practice, the measured signals are limited in time and the FFT calculates the frequency transform over a certain number of discreet frequencies called bins.

Spectral Leakage:

In reality, signals are of time-limited nature and nothing can be known about the signal beyond the measured interval. For example, if the measurement of a never ending continuous train of sinusoidal wave is of interest, at some point of time we need to terminate our measurement to do further analysis. The limit on the time is also posed by limitations of the measurement system itself (like buffer size) besides other factors. Some assumptions have to made about the nature of the signal outside the measured interval. Fourier Transforms implicitly assumes that the signal essentially repeats itself after the measured interval.

Figure 1 illustrates the scenario in which a continuous train of sinusoidal signal is observed over a finite interval of time ("measured signal"). As discussed, the FFT assumes the signal to be continuous (conceptually, it does this by juxtaposing the measured signal repetitively). Observe the glitches in the assumed signal. These glitches are the manifestations of the measurement time relative to the frequency of the actual signal. If measurement time is an integral multiple of the rate of the actual signal (i.e. the inverse of the frequency of the signal), then no glitch will be observed in the assumed signal. In Figure 1, the measurement time is purposefully made to be a non-integral multiple of the actual signal rate. These sharp discontinuities will spread out in the frequency domain. This is called spectral leakage.

Figure 1: Impact of observation interval on FFT

Lets visualize the concept of spectral leakage by taking a pure sinusoidal signal as an example. Here we consider a sinusoid of 7Hz frequency (7 cycles in 1 second) and sampling it with a sampling frequency Fs=100Hz. We observe the signal for 100 seconds (700 cycles in total) and take the FFT of the observed signal. Figure 2 illustrates the frequency spectrum of the observed/measured signal. Essentially, the frequency spectrum contains a distinct peak at 7Hz.

Figure 2: Frequency Spectrum of a 7 Hz sinusoid observed for 100 seconds

Next, the measurement window is shrunk to 1 second. Now the FFT is taken on this observed signal. Figure 3 illustrates the frequency spectrum of the observed signal (which has 7 cycles only). The frequency spectrum contains some spectral leakage because of limited observation interval.

Figure 3 : A 7 Hz sinusoid observed for 1 second and its Frequency Spectrum

In the previous illustrations (Figure 2 and 3), the observation time interval contained an integral number of sinewave cycles (i.e. in Figure 2 the observation interval contained exactly 700 cycles of sinusoid and in Figure 3 the observation interval contained exactly 7 cycles of sinusoid in the time domain).

For the next illustration the measurement time interval is adjusted in such a way that number of cycles in the observation window is no longer an integer. In Figure 4, the signal is observed for 1.4 seconds, which implies that there are 9.8 cycles. Now the observed signal does not end at zero amplitude at t=1.4 seconds. This scenario gives rise to gliches in the FFT's assumed signal (which constructs a periodic signal from the observed signal) and obviously results in more spectral leakage (Compare the Frequency spectrum in Figure 3 and Figure 4).

Figure 4 : A 7 Hz sinusoid observed for 1.4 seconds and its Frequency Spectrum.

The effect of spectral leakage may be lessened if the observed signal does not contain any discontinuity at the end of the measurement time (this scenario rarely occurs in any real application). Another scenario in which the spectral leakage can be reduced is by having a signal that gradually reduces to zero at the ends of the measurement time.( All the windows (like Hamming, Hanning, Bartlett, etc.., essentially attempts to do this). Such signal would have no discontinuity when made periodic and so does not suffer spectral leakages.

Conclusion:

Spectral leakage is not due to FFT but due to the finite observation time. Spectral Leakage gives rise to two problems : 1) The spectral component of the desired signal no longer contains the complete energy, rather it also contains the energy of adjacent components and noise, thereby reducing the Signal to noise ratio; 2) Spectral leakage from a larger signal component may significantly overshadow other smaller signals making them difficult to identify or detect.

Windowing Techniques are used to mitigate the effects of spectral leakage and therefore the restriction of having a limited observation interval.

Matlab Code:

Here is a matlab code to simulate the charts given in this discussion

%Illustration of Spectral Leakage and its effects on Frequency Spectrum

%(FFT)

%Author: Mathuranathan Viswanathan

%http:gaussianwaves.blosgpot.com

%License - Creative Commons - cc-by-nc-sa

observationTime=1; %Input : Observation time interval change it to 100 and 1.4

%and see the effect on FFT

Fx=7; %Frequency of the sinusoid

Fs=100; %Sampling Frequency
t=0:1/Fs:observationTime;
x=1*sin(2*pi*Fx*t);
plot(t,x)
title('SineWave - Frequency = 7Hz');
xlabel('Times (s)');ylabel('Amplitude');

%Perform FFT

NFFX=2.^(ceil(log(length(x))/log(2)));
FFTX=fft(x,NFFX);%pad with zeros
NumUniquePts=ceil((NFFX+1)/2);
FFTX=FFTX(1:NumUniquePts);
MY=abs(FFTX);
MY=MY*2;
MY(1)=MY(1)/2;
MY(length(MY))=MY(length(MY))/2;
MY=MY/length(x);
f1=(0:NumUniquePts-1)*Fs/NFFX;

figure;
stem(f1,MY,'k');
title('FFT of the sine wave');
axis([0 50 0 1]);
xlabel('Frequency');ylabel('Amplitude');

Reference:

[1] Fredric J Harris, "On the use of windows for Harmonic Analysis with the Discrete Wavelet Transform", Proceedings of IEEE, Vol 66,No 1, January 1978 - click here

Recommended Readings:

Raised Cosine Filter

2010-12-27T04:09:00.001Z

Role of Pulse Shaping in Communication Systems:
In communication systems, data is transmitted as binary bits ( ones and zeros). It is easier to implement a binary system using switches, where turning on a switch represents '1' and turning it off represents '0'. Such simple binary systems essentially represents ones and zeros as rectangular pulses of finite duration (say τ).
In frequency domain a rectangular pulse of finite duration τ manifests itself as a sinc pulse of infinite duration. Also, most of the energy of the rectangular pulse is concentrated with in -1/τ to 1/τ in the frequency domain. This implies a pulse of duration τ requires twice its bandwidth for reliable transmission. This poses a limitation when it comes to bandwidth limited system, if we wish to increase the transmission data rate. The challenge here is, how do we transmit the data with highest possible data rate in a given band-limited system. The answer is : Pulse Shaping.

In a band-limited system, when we try to increase the data rate, it may lead to Inter Symbol Interference (ISI). There are two criteria for a non-interference systems where pulse shaping is employed. (1) The pulse shape exhibits a zero crossing at the sampling point of all pulse intervals except its own ; (2) the shape of the pulses be such that the amplitude decays rapidly outside of the pulse interval. A rectangular pulse satisfies the first criterion (where it contains zero crossing - see figure below) but not the second criterion ( the energy of the rectangular pulse does not decay rapidly outside the pulse interval and in fact it extends to infinite bandwidth). A raised cosine pulse is designed to satisfy these two criterion and thus provides a ISI free system.

Raised Cosine Filters/Pulses:
A Raised Cosine pulse takes on the shape of a sinc pulse in time domain and is given by the following function (This equation is apt for digital domain and matlab simulation, it is obtained from its analog form by substituting "t" by n*Ts)

Here Ts is the Sampling period ,n is the sample number, α is a parameter that governs the bandwidth occupied by the pulse and the rate at which the tails of the pulse decay. A value of α = 0 offers the narrowest bandwidth, but the slowest rate of decay in the time domain. When α = 1,the bandwidth is 1/τ, but the time domain tails decay rapidly. For more details on investigation of raised cosine characteristics and its implementation - refer [1]

The following Matlab program simulated a Raised Cosine Filter and plots its frequency response and its eye diagram. (Click here to download)

L=41; %filter Length
R=1E6; %data Rate = 1Mbps
Fs=8*R; %oversampling by 8
T=1/R; %pulse duration
Ts=1/Fs; %sampling duration
alpha =0.5; % Design Factor for Raised Cosing Filter
%----------------------------------------------------------
%Raised Cosing Filter Design
%----------------------------------------------------------
if mod(L,2)==0
        M=L/2 ; % for even value of L
else
         M=(L-1)/2; % for odd value of L
end
g=zeros(1,L); %Place holder for RC filter's transfer function
for n=-M:M
        num=sin(pi*n*Ts/T)*cos(alpha*pi*n*Ts/T);
       den=(pi*n*Ts/T)*(1-(2*alpha*n*Ts/T)^2);
        g(n+M+1)=num/den;
        if (1-(2*alpha*n*Ts/T)^2)==0
                g(n+M+1)=pi/4*sin(pi*n*Ts/T)/(pi*n*Ts/T);
        end
        if (pi*n*Ts/T)==0
                g(n+M+1)=cos(alpha*pi*n*Ts/T)/(1-(2*alpha*n*Ts/T)^2);
        end
end
%----------------------------------------------------------
% Plot the transfer function of RC filter
figure;
impz(g,1);

The above section of code implements the raised cosine filter equation and plots the impulse response of the filter.The RC filtered outputs for α = 1 and α =0.5 are plotted below. As we can see that as α increases (from 0 to 1), the pulse duration decreases (in frequency domain it will consume excess bandwidth). We must strict an optimum balance between the value of α and the bandwidth it occupies, in order for the pulse to travel undistorted in a given band-limited channel.

Next we generate a random data stream and upsample it with an upsample factor of Fs/R (where Fs is the sampling rate and R is the data rate). The upsampled data is convolved with the designed filter to give the RC filtered output.

%----------------------------------------------------------

%Generate data of random 1s and 0s

data=2*(rand(1,1000)>=0.5)-1; %Polar encoding -> 1-> 1V, 0->-1V
output=upsample(data,Fs/R);
%y=conv(g,output); %Convolving the data signal with the Raised Cosine Filter
y=filter(g,1,output); %you can use either Conv function or filter function to obtain the output

%----------------------------------------------------------

%Plot data and RC filtered Output

figure;
subplot(2,1,1);
stem(data);
title('Input data to the Raised Cosine Filter');
xlabel('Samples');
ylabel('Amplitude');
subplot(2,1,2);
plot(y);
title('Response of the Raised Cosine Filter for the given Input');
xlabel('Samples');
ylabel('Amplitude');

The next section of code computes the frequency response (by employing FFT) and also plots the eye diagram for the RC filtered output signal.

%Obtain FFT of the output to plot its frequency response.

Fn=Fs/2;NFFY=2.^(ceil(log(length(y))/log(2)));
FFTY=fft(y,NFFY);%pad with zeros
NumUniquePts=ceil((NFFY+1)/2);
FFTY=FFTY(1:NumUniquePts);
MY=abs(FFTY);
MY=MY*2;
MY(1)=MY(1)/2;
MY(length(MY))=MY(length(MY))/2;
MY=MY/length(y);
f1=(0:NumUniquePts-1)*2*Fn/NFFY;

%Plot Frequency spectrum

figure;
plot(f1,20*log10(abs(MY).^2));xlabel('FREQUENCY(Hz)');ylabel('DB');
grid
title('Frequency domain plots')

% Plot Eye Diagram

plotEye(y,Fs/R);

The frequency response and the eye-diagram for the Raised Cosine Filtered output for α=1 and α=0.5 are plotted below for comparison. The frequency response shows minimized side lobes for α=1, whereas for the case α=0.5 the side lobes' energy is higher when compared to α=1. Also, for α=1, the eye opening is clearer when compared to α=0.5. This implies that as we increase the value of α from 0 to 1, the ISI decreases where as the excess bandwidth of the signal increases. An optimum balance between the value of α and the bandwidth of the band-limited channel is necessary for ISI free reliable communication between the transmitter and receiver.

Note : The function plotEye is a custom made function which does not depend on the Matlab Communication toolbox. Please click here to obtain the code for this function.

Reference:
Ken Gentile,"The care and feeding of digital, pulse-shaping filters",RF design.com magazine,April 2002.

Recommended Products/books on Matlab:

An utility for plotting eyediagram in matlab (without using Communication toolbox)

2010-12-26T06:43:00.006Z

Here is a small function written in Matlab to plot eyediagram for any real signal. This function will come handy when the communication toolbox of Matlab is not accessible/installed. This function will work with basic installation of Matlab.

Syntax:

plotEye(input,overSamplingRate)
where,
input = input signal
overSamplingRate = oversampling factor used to sample the signal. Generally this must be Fs/R where Fs is the sampling rate and R is the data rate.

We will use this routine in subsequent post to analyze Inter Symbol Interference (ISI)

Matlab Code: (click here to download)

function [T,EyeSig]=plotEye(inputSignal,overSampleRate,offset);
%----------------------------------------------
%Function to plot eye diagram of inputSignal
%Inputs
% inputSignal-real signal
% overSampleRate-oversample rate
% http://gaussianwaves.blogspot.com
% License: Creative Commons - CC - BY-NC-SA
%----------------------------------------------
if(nargin < 3),offset=0;end
inputSignal=interp(inputSignal,4);
overSampleRate=overSampleRate*4;
f=mod(length(inputSignal),overSampleRate);
inputSignal=real(inputSignal(f+1+offset:end-(overSampleRate-offset)));
%------------------------------------------------
EyeSigRef=reshape(real(inputSignal),overSampleRate,length(inputSignal)/overSampleRate);
EyeSigm1=EyeSigRef(overSampleRate/2+2:overSampleRate,1:end-2);
EyeSig0=EyeSigRef(:,2:end-1);
EyeSigp1=EyeSigRef(1:overSampleRate/2,3:end);
EyeSig=[EyeSigm1;EyeSig0;EyeSigp1];

T=[-1+1/(overSampleRate):1/(overSampleRate):+1-1/(overSampleRate)];
figure;
plot(T,EyeSig),xlabel('Symbol periods'),ylabel('Amplitude');

Sample Plot:

Recommended books/products on Matlab:

Moving Average Filter ( MA filter )

2010-11-23T07:15:00.004Z

The moving average filter is a simple Low Pass FIR (Finite Impulse Response) filter commonly used for smoothing an array of sampled data/signal. It takes M samples of input at a time and take the average of those M-samples and produces a single output point. It is a very simple LPF structure that comes handy for scientists and engineers to filter unwanted noisy component from the intended data.

As the filter length increases ( the parameter M ) the smoothness of the output increases, whereas the sharp transitions in the data are made increasingly blunt. This implies that this filter has excellent time domain response but a poor frequency response.

The MA filter perform three important functions:
1) It takes M input points, computes the average of those M-points and produces a single output point
2) Due to the computation/calculations involved , the filter introduces a definite amount of delay
3) The filter acts as a Low Pass Filter (with poor frequency domain response and a good time domain response).

Matlab Code: (click here to download)
Following matlab code simulates the time domain response of a M-point Moving Average filter and also plots the frequency response for various filter lengths.

%Moving Average Filter(MA) Demonstration
%Computes and plots the time domain output of a M-point Moving Average
%Filter. Note : The parameter M must be a odd number inorder to realize the
%filter.
%Plots the frequency response of the M-point Moving Average filter for
%comparison.
%Author : Mathuranathan Viswanathan
%http://gaussianwaves.blogspot.com
%License - Creative Commons - cc-by-nc-sa

clear;
clc;

M=51; %Number of taps/points in the MA filter
data=[zeros(1,1000),5*ones(1,3000),zeros(1,1000)]; %5000 points sample data
noise= 2*rand(1,5000)-1;

x=data+noise; %Adding some noise to the data

%Plot the input
figure;
plot(x);
xlabel('Samples');
ylabel('Amplitude');
title('Input to the M-point Moving Average Filter');

%Computing the output response of the given moving average filter
figure;
y=zeros(1,5000);
for i=(M-1)/2+1:length(data)-(M-1)/2
    y(i)=0;
    for k=-(M-1)/2:(M-1)/2
        y(i)=y(i)+x(i+k);
    end
    y(i)=y(i)/M;
end
plot(y);
xlabel('Samples');
ylabel('Amplitude');
title(['Filtered output of the ',num2str(M),'-point Moving Average filter']);

%========================================================================
%Frequency Response for various Filter lengths
%========================================================================
colors=['r','b','k','g']; %To plot in different colors
colorIndex=1;
figure;
M=[3,11,51,101];
for i=M
    b=repmat(1/i,1,i);
    a=1; %Since this is a FIR filter a=1 and b=[1/M,1/M,1/M,...,repeated M times]
    [h,w]=freqz(b,a); %Using the Freqz function to get the frequency response
    h=abs(h);
    w=w/max(w);
    plot(w,h,colors(colorIndex)); %Plotting the response
    colorIndex=colorIndex+1;
    hold on;
end
hold off;
legend('M=3','M=11','M=51','M=101');
xlabel('Normalized Frequency');
ylabel('Magnitude');
title('Frequency Response of M-point Moving average filters');

Time domain Response:

Input to MA filter	3-point MA filter output

51-point MA filter output	101-point MA filter output

501-point MA filter output

On the first plot, we have the input that is going into the moving average filter. The input is noisy and our objective is to reduce the noise. The next figure is the output response of a 3-point Moving Average filter. It can be deduced from the figure that the 3-point Moving Average filter has not done much in filtering out the noise. We increase the filter taps to 51-points and we can see that the noise in the output has reduced a lot, which is depicted in next figure. We increase the taps further to 101 and 501 and we can observe that even-though the noise is almost zero, the transitions are blunted out drastically (observe the slope on the either side of the signal and compare them with the ideal brick wall transition in our input).

Frequency Response:

From the frequency response it can be asserted that the roll-off is very slow and the stopband attenuation is not good. Given this stop band attenuation,clearly, the moving average filter cannot separate one band of frequencies from another. As we know that a good performance in the time domain results in poor performance in the frequency domain, and vice versa. In short, the moving average is an exceptionally good smoothing filter (the action in the time domain), but an exceptionally bad low-pass filter (the action in the frequency domain).

For Further Reading:

Random Interleaver

2010-11-16T09:15:00.002Z

Random Interleaver:

The Random Interleaver rearranges the elements of its input vector using a random permutation. The incoming data is rearranged using a series of generated permuter indices. A permuter is essentially a device that generates pseudo-random permutation of given memory addresses. The data is arranged according to the pseudo-random order of memory addresses.

The de-interleaver must know the permuter-indices exactly in the same order as that of the interleaver. The de-interleaver arranges the interleaved data back to the original state by knowing the permuter-indices.

The interleaver depth (D) (Not sure what this term means ? - check out this post - click here) is essentially the number of memory addresses taken for permutation at a time. If the number of memory addresses taken for permutation increases, the interleaver depth increases.

In the Matlab simulation that is given below, an interleaver depth of 15 is used for illustration. This means that 15 letters are taken at a time and are permuted (rearranged randomly). This process is repeated consecutively for the next block of 15 lettes.

As you may observe from the simulation that increasing the interleaver depth will increase the degree of randomness in the interleaved data and will decrease the maximum burst length after the de-interleaver operation

Matlab Code: (click here to download)

A sample matlab code that simulates the above mentioned random interleaver is given below. The input data is a repeatitive stream of following symbols - "THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_". This code simulates only the interleaving/denterleaving part.The burst errors produced by channel are denoted by '*'.
.

%Demonstration of PseudoRandom Interleaver/Deinterleaver
%Author : Mathuranathan for http://gaussianwaves.blogspot.com
%License - Creative Commons - cc-by-nc-sa

clc;
clear;
%____________________________________________
%Input Parameters
%____________________________________________
D=15; %Interleaver Depth. More the value of D more is the randomness
n=5; %Number of Data blocks that you want to send.
data='THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_'; % A constant pattern used as a data
%____________________________________________

%Make the length of the data to be a multiple of D. This is for
%demonstration only.Otherwise the empty spaces has to be zero filled.
data=repmat(data,1,n); %send n blocks of specified data pattern
if mod(length(data),D) ~=0
    data=[data data(1:(D*(fix(length(data)/D+1))-length(data)))];
end

%We sending D blocks of similar data
%intlvrInput=repmat(data(1:n),[1 D]);

%fprintf('Input Data to the Interleaver -> \n');
%disp(char(intlvrInput));
%disp('____________________________________________________________________________');

%INTERLEAVER
%Writing into the interleaver row-by-row
permuterIndex=randperm(D);
intrlvrOutput=[];
index=1;
for i=1:fix(length(data)/D)
    intrlvrOutput=[intrlvrOutput data(permuterIndex+(i-1)*D)];
end
for i=1:mod(length(data),D)
    intrlvrOutput=[intrlvrOutput data(permuterIndex(i)+fix(length(data)/D)*D)];
end

uncorruptedIntrlvrOutput=intrlvrOutput;

%Corrupting the interleaver output by inserting 10 '*'
intrlvrOutput(1,25:34)=zeros(1,10)+42;

%DEINTERLEAVER
deintrlvrOutput=[];
for i=1:fix(length(intrlvrOutput)/D)
    deintrlvrOutput(permuterIndex+(i-1)*D)=intrlvrOutput((i-1)*D+1:i*D);
end
for i=1:mod(length(intrlvrOutput),D)
    deintrlvrOutput((fix(length(intrlvrOutput)/D))*permuterIndex+i)=intrlvrOutput((i+1):(i+1)*D);
end
deintrlvrOutput=char(deintrlvrOutput);

disp('Given Data -->');
disp(data);
disp(' ')

disp('Permuter Index-->')
disp(permuterIndex);
disp(' ')

disp('PseudoRandom Interleaver Output -->');
disp(uncorruptedIntrlvrOutput);
disp(' ')

disp('Interleaver Output after being corrupted by 10 symbols of burst error - marked by ''*''->');
disp(char(intrlvrOutput));
disp(' ')

disp('PseudoRandom Deinterleaver Output -->');
disp(deintrlvrOutput);

Output

Given Data -->
THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_Q

Permuter Index-->
15 1 11 9 7 12 5 6 13 4 2 10 3 8 14

PseudoRandom Interleaver Output -->
NTBKIRQUO_H_ECWR__PUO_JVXFSOMET_DYAO_LGET_HZ__HR_COUIWQEB_KN_FOSMVJUE_O_XPRHTO_ZGLA__HDEYTFEOBKWICNU_RQ__TOV_PEUMRJXO_S_EHGDY_AZTLEO__HO_WR_NCK_IQOUBFHXEOSRMP_U_VJ_T_E_O_TZYHA_GLDEXQNOB_K_FCUWIROE_RV__PSTMJEUOHQ_TGDHY_EZL_AO_

Interleaver Output after being corrupted by 10 symbols of burst error - marked by '*'->
NTBKIRQUO_H_ECWR__PUO_JV**********AO_LGET_HZ__HR_COUIWQEB_KN_FOSMVJUE_O_XPRHTO_ZGLA__HDEYTFEOBKWICNU_RQ__TOV_PEUMRJXO_S_EHGDY_AZTLEO__HO_WR_NCK_IQOUBFHXEOSRMP_U_VJ_T_E_O_TZYHA_GLDEXQNOB_K_FCUWIROE_RV__PSTMJEUOHQ_TGDHY_EZL_AO_

PseudoRandom Deinterleaver Output -->
THE_QUICK_BROWN_***_JU*P*_OV*R*THE_LAZ*_*OG_*HE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_Q

As we can see from the above simulation that, eventhough the channel introduces 10 symbols of consecutive burst error, the interleaver/deinterleaver operation has effectively distributed the errors and reduced the maximum burst length to 3 symbols.

Recommended Books

Block Interleaver Design for RS codes

2010-10-28T09:18:00.001+01:00

Introduction:
A (n,k) Reed Solomon (RS) encoder, takes k user data symbols and converts it into a n symbol wide codeword, by adding n-k parity symbols. The error correcting capability (t) of the RS code is computed as t<=(n-k)/2. That is, a RS code with n-k parity symbols can correct a burst error of upto (n-k)/2 symbol errors.

Block Interleavers:
Suppose, assume that the dominant error mechanism in a channel is of burst type. A burst of length b is defined as a string of b unreliable consecutive symbols. If the expected burst length, b is less than or equal to t (the number of correctable symbol errors by RS coding), the code can be used as it is. However, if bursts length b >t, the error correcting code will fail. This is where interleaving comes to our rescue.

Let us assume that b > t and b ≤ t × d, where d (the interleaving depth) is an integer. The Reed-Solomon (n,k) code can be used if we can spread the burst error sequence over several code blocks so that each block has no more than t errors (which can then be corrected). This can be accomplished using block interleaving as follows. Instead of encoding blocks of k symbols and then sending the encoded symbols consecutively, we can interleave the encoded blocks and transmit the interleaved data. In the case where n = 255, k = 235, t = 10, d = 5,the data bytes output from the Reed-Solomon encoder would appear as shown below , where bytes numbered 0 to 234 are the data bytes and bytes 235 to 254 are the parity check bytes.

Here, the data is written row by row and read back column by column.Consider now the effect of a burst error of length b > t, (where t is the number of correctable errors per block) and b ≤ v × d for some v, on the received symbols in the table. Because of the order in which the symbols are sent, a burst length less than or equal to d × i will effect at most d + 1 consecutive columns of the table, depending on where the burst
starts. Notice that any single row (which corresponds to a codeword) has no more than v errors in it. If d > t, these errors are within the error correction capability of the code and can be corrected. In this case, d becomes the interleaving depth. The trade-off is that extra buffer space is required to store the interleaver table and additional delay is introduced. The worst case burst length determines the size of the table (and the interleaving depth) and the table size determines the amount of buffer space required and the delay.

Design Example:

Consider a (255,235) Reed Solomon coding system. This code can correct upto t=(n-k)/2 = 10 symbols errors. Lets assume that the channel that we are going to use, is expected to cause b=253 symbols. Then the interleaver depth (d) is calculated as

In this case , an interleaver depth of 26 is enough to combat the burst errors introduced by the channel. The block interleaver dimensions would be d x n = 26 x 255 (26 rows by 255 columns).

Matlab Code: (click here to download)

A sample matlab code that simulates the above mentioned block interleaver design is given below. The input data is a repeatitive stream of following symbols - "THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_". A (255,235) Reed Solomon decoder (with correction capability of 10 symbols) is used. We assume that the channel is expected to produce a maximum of consecutive 20 symbols of burst error. The burst errors are denoted by '*'.
.

%Demonstration of design of Block Interleaver for Reed Solomon Code
%Author : Mathuranathan for http://gaussianwaves.blogspot.com
%License - Creative Commons - cc-by-nc-sa
clc;
clear;
%____________________________________________
%Input Parameters
%____________________________________________
%Interleaver Design for Reed-Solomon Codes
%RS code parameters
n=255; %RS codeword length
k=235; %Number of data symbols
b=20; %Number of symbols that is expected to be corrupted by the channel
%____________________________________________

p=n-k; %Number of parity symbols
t=p/2; %Error correction capability of RS code

fprintf('Given (%d,%d) Reed Solomon code can correct : %d symbols \n',n,k,fix(t));
fprintf('Given - expected burst error length from the channel : %d symbols \n',b);
disp('____________________________________________________________________________');
if(b>t)
    fprintf('Interleaving MAY help in this scenario\n');
else
    fprintf('Interleaving will NOT help in this scenario\n');
end
disp('____________________________________________________________________________');
D=ceil(b/t)+1; %Intelever Depth

memory = zeros(D,n); %constructing block interleaver memory

data='THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_'; % A constant pattern used as a data

%If n>length(data) then repeat the pattern to construct a data of length 'n'
data=char([repmat(data,[1,fix(n/length(data))]),data(1:mod(n,length(data)))]);

%We sending D blocks of similar data
intlvrInput=repmat(data(1:n),[1 D]);

fprintf('Input Data to the Interleaver -> \n');
disp(char(intlvrInput));
disp('____________________________________________________________________________');

%INTERLEAVER
%Writing into the interleaver row-by-row
for index=1:D
    memory(index,1:end)=intlvrInput((index-1)*n+1:index*n);
end
intlvrOutput=zeros(1,D*n);
%Reading from the interleaver column-by-column
for index=1:n
    intlvrOutput((index-1)*D+1:index*D)=memory(:,index);
end

%Create b symbols error at 25th Symbol location for test in the interleaved output
%'*' means error in this case
intlvrOutput(1,25:24+b)=zeros(1,b)+42;
fprintf('\nInterleaver Output after being corrupted by %d symbol burst error - marked by ''*''->\n',b);
disp(char(intlvrOutput));
disp('____________________________________________________________________________');

%Deinteleaver
deintlvrOutput=zeros(1,D*n);
%Writing into the deinterleaver column-by-column
for index=1:n
    memory(:,index)=intlvrOutput((index-1)*D+1:index*D)';
end

%Reading from the deinterleaver row-by-row
for index=1:D
    deintlvrOnput((index-1)*n+1:index*n)=memory(index,1:end);
end
fprintf('Deinterleaver Output->\n');
disp(char(deintlvrOnput));
disp('____________________________________________________________________________');

Output

Given (255,235) Reed Solomon code can correct : 10 symbols
Given - expected burst error length from the channel : 20 symbols
____________________________________________________________________________
Interleaving MAY help in this scenario
____________________________________________________________________________
Input Data to the Interleaver ->
THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_____________________________________________________________________________

Interleaver Output after being corrupted by 20 symbols burst error - marked by '*'->
TTTHHHEEE___QQQUUUIIICCC********************N___FFFOOOXXX___JJJUUUMMMPPPSSS___OOOVVVEEERRR___TTTHHHEEE___LLLAAAZZZYYY___DDDOOOGGG___TTTHHHEEE___QQQUUUIIICCCKKK___BBBRRROOOWWWNNN___FFFOOOXXX___JJJUUUMMMPPPSSS___OOOVVVEEERRR___TTTHHHEEE___LLLAAAZZZYYY___DDDOOOGGG___TTTHHHEEE___QQQUUUIIICCCKKK___BBBRRROOOWWWNNN___FFFOOOXXX___JJJUUUMMMPPPSSS___OOOVVVEEERRR___TTTHHHEEE___LLLAAAZZZYYY___DDDOOOGGG___TTTHHHEEE___QQQUUUIIICCCKKK___BBBRRROOOWWWNNN___FFFOOOXXX___JJJUUUMMMPPPSSS___OOOVVVEEERRR___TTTHHHEEE___LLLAAAZZZYYY___DDDOOOGGG___TTTHHHEEE___QQQUUUIIICCCKKK___BBBRRROOOWWWNNN___FFFOOOXXX___JJJUUUMMMPPPSSS___OOOVVVEEERRR___TTTHHHEEE___LLLAAAZZZYYY___DDDOOOGGG___TTTHHHEEE___QQQUUUIIICCCKKK___BBBRRROOOWWWNNN___FFFOOOXXX___JJJUUUMMMPPPSSS___OOOVVVEEERRR___TTTHHHEEE_______________________________________________________________________________
Deinterleaver Output->
THE_QUIC*******_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_THE_QUIC*******_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_THE_QUIC******N_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_LAZY_DOG_THE_QUICK_BROWN_FOX_JUMPS_OVER_THE_____________________________________________________________________________

As we can see from the above simulation that, eventhough the channel introduces 20 symbols of consecutive burst error (which is beyond the correction capability of the RS decoder), the interleaver/deinterleaver operation has effectively distributed the errors and reduced the maximum burst length to 7 symbols (which is easier to correct by (255,235) Reed Solomon code.

Recommended Books:

Interleavers and deinterleavers

2010-10-25T11:36:00.004+01:00

Interleavers and Deinterleavers are designed and used in the context of characteristics of the errors that might occur when the message bits are transmitted through a noisy channel. To understand the functions of an interleaver/deinterleaver, understanding of error characteristics is essential. Two types are errors concern communication system design engineer. They are burst error and random error

Random Errors:

Error locations are independent of each other. Error on one location will not affect the errors on other locations. Channels that introduce these types of errors are called channels without memory (since the channel has no knowledge of error locations since the error on location does not affect the error on another location).

Burst Errors:

Errors are depended on each other. For example, in channels with deep fading characteristics, errors often occur in bursts (affecting consecutive bits). That is, error in one location has a contagious effect on other bits. In general, these errors are considered to be dependent and such channels are considered to be channels with memory.

Interleaver/Deinterleaver :

One of the most popular ways to correct burst errors is to take a code that works well on random errors and interleave the bursts to "spread out" the errors so that they appear random to the decoder. There are two types of interleavers commonly in use today, block interleavers and convolutional interleavers.

The block interleaver is loaded row by row with L codewords, each of length n bits. These L codewords are then transmitted column by column until the interleaver is emptied. Then the interleaver is loaded again and the cycle repeats. At the receiver, the codewords are deinterleaved before they are decoded. A burst of length L bits or less will cause no more than 1 bit error in any one codeword. The random error decoder is much more likely to correct this single error than the entire burst.The parameter L is called the interleaver degree, or interleaver depth. The interleaver depth is chosen based on worst case channel conditions. It must be large enough so that the interleaved code can handle the longest error bursts expected on the channel. The main drawback of block interleavers is the delay introduced with each row-by-row fill of the interleaver.

Convolutional interleavers eliminate the problem except for the delay associated with the initial fill. Convolutional interleavers also reduce memory requirements over block interleavers by about one-half [1]. The big disadvantage of either type of interleaver is the interleaver delay introduced by this initial fill. The delay is a function of the interleaver depth and the data rate and for some channels it can be several seconds long. This long delay may be unacceptable for some applications. On voice circuits (as in GSM), for example, interleaver delays confuse the unfamiliar listener by introducing long pauses between speaker transitions. Even short delays of less than one second are sufficient to disrupt normal conversation. Another disadvantage of interleavers is that a smart jammer can choose the appropriate time to jam to cause maximum damage. This problem is overcome by randomizing the order in which the interleaver is emptied.

In practice, interleaving is one of the best burst-error correcting techniques. In theory, it is the worst way to handle burst errors. Why? From a strict probabilistic sense, we are converting "good" errors into "bad" errors. Burst errors have structure and that structure can be exploited. Interleavers "randomize" the errors and destroy the structure. Theory differs from reality,however. Interleaving may be the only technique available to handle burst errors successfully.

For example, Viterbi [2] shows that, for a channel impaired by a pulse jammer, exploiting the burst structure is not enough. Interleaving is still required. This does not mean that we should be careless about our choice of code and take up the slack with long interleavers. Codes designed to correct burst errors can achieve the same performance with much shorter interleavers. Until the coding theorists discover a better way, interleaving will be an essential error control coding technique for bursty channels.

References :

[1]B. Sklar, Digital Communications Fundamentals and Applications, Englewood Cliffs, New Jersey: Prentice Hall, 1988.
[2]A. J. Viterbi, "Spread Spectrum Communications – Myths and Realities," IEEE Communications Magazine, vol. 17, No. 5, pp. 11-18, May 1979.

Recommended Books: