Array pattern multiplication – phased array antenna

Articles in this series

Array pattern multiplication

If the element pattern of an array is very broad (as in an Omni directional antenna element) the pattern of the array is mostly the Array Factor since it is more directional. However, if the element pattern has some distinctive features by employing specific illumination function, then the composite radiation pattern is obtained by the principal of pattern multiplication. If Farrayis the Fourier transform of the array field and Felement is the Fourier transform of the element field, the combined pattern is given by the pattern multiplication

Composite antenna pattern = Farray × Felement

Example:

  1. Uniformly illuminated array produces a first sidelobe level of -13.5 dB [=20×log10 (field amplitude)]. Due to the vast difference in levels between major and minor lobes, logarithmic (dB) scale is preferred in such measurements to provide the required dynamic range. This illumination function produces a better aperture efficiency (gain), sharper beamwidth but at the cost of near sidelobes getting elevated.
Radiation pattern with uniform aperture illumination
Figure 1: Radiation pattern with uniform aperture illumination

2. The cosine aperture distribution on a pedestal attempts to deal with a maximized gain along with control on the sidelobe levels. The accompanying figure below is the array pattern with illumination function of the type f(x) =0.9×cos(x) +0.1, where x refers to the linear position of the element.

Figure 2: Radiation pattern with cosine aperture illumination

Instead of tapering the edge illumination to zero, the cosine function sits on a pedestal to provide a constant illumination, but at a lower percentage level, at the peripheral geometry of the radiating element. The main beam is broadened but the first sidelobes are at -23.5dB level, an order of improvement on the -13.5 dB level of the previous case.

Several sophisticated illumination functions are available to the designer now to optimize array design.

Concluding remarks

It should be stressed that a Phased Array produces huge opportunity for designers to exploit its dynamic properties; but at the same time, it challenges us with system engineering and integration complexity, maintainability and finally in the cost factor. In an evolving era of technology, limitations in engineering and operation, would always witness a better tomorrow, provided the basic scientific validation and sustainability are present in all of the newer innovations. 

Rate this post: Note: There is a rating embedded within this post, please visit this post to rate it.

Grating Lobes in Electronic Scanning

Note: There is a rating embedded within this post, please visit this post to rate it.

The array is a periodic structure in space with elements placed at equidistant. d/λ acts as the frequency term and λ/d as the time period. The incremental phase shift between adjacent elements is (2π/λ)×d×[sin(θ)-sin(θ0)]d is the inter-element distance in unit of λ, θ-any position during the scanning of the array and θ0 – the angular direction set for the array to transmit/receive with the main lobe having the maximum gain.

For the ensuing discussion here, let us assume that no preferred direction θ0 is set and the array is free to scan from its broadside (θ=0°) to the limits of scan ± 90° or u=± 1.  The inter-element phase shift now is (2π/λ) × d × u, where u = sin θ, and the magnitude of phase shift attains a maximum at |u|=1.

The below Fig.1 traces the phase function for various values of d, as it is increased from 0.25 λ to 2 λ.

Fig. 1 Phase range Vs. d/λ

Two apparent factors emerge from this plot: As d is increased between the elements, the rate of phase change is higher and the range between u=-1 and u=1 is also increased. The field is sampled (for its amplitude and phase) at each element and the sampling is repeated at the spatial interval of d. Up to a distance of d≤ λ/2, the maximum phase difference between adjacent elements is ≤ π for |u|=1, and increases beyond this limit for d> λ/2. This means that samples between elements are taken at least once, within a phase difference of π, as long as d≤ λ/2.  The visible phase region for a scan limit from u=-1 to u= 1, is when d≤ λ/2 condition exists; and taking at least one sample per π radians satisfies the Nyquist criterion (π radians being the Nyquist space). This condition no longer exists when the inter-element distance is increased beyond λ/2.

For example, if d= λ, the maximum phase difference between adjacent elements is for |u|=1. But the sampling once between adjacent elements now spans a phase range of , resulting in a clear case of under sampling in the spatial domain. The situation becomes even more critical when d is further increased. As observed earlier, the rate of phase difference between elements and the total range for (u=-1 to u=1), climb up to such levels, much poorer sampling rates happen among the elements. This results in increasing spatial aliasing, akin to what one experiences in the case of under sampled signals in the frequency domain.

Such spatial aliasing result in the occurrence of extra peaks (with near similar values of the main beam) and thus create disruptions in the antenna pattern. They result in giving wrong directional readings as well as sharing the power and distributing it in unwanted zones of the antenna array pattern. Such events are named as the presence of grating lobes.

Examination of the occurrence and placement of grating lobes:

 As seen in the earlier section, the resultant field of a linear array of N elements is given by:

[ E= the field intensity at each element (transmit/receive) and δ=(2π/λ)×d×(u-u0) ]

Grating lobes were discussed so far only with respect to the increase of inter-element distance > λ/2. But there are other contributing factors too.

  1. The inter-element phase δ has three variables such as λ, d and u0. For narrow band systems, the dispersion in the operating wavelength is small and may be ignored; but in wideband and ultra-wideband systems it is a factor to reckon with. In such cases, both phase shift and time delay circuits will have to be employed.
  2. The effect of increasing d >λ/2 has already been discussed for the occurrence of spatial aliasing.
  3. The setting of u0 has limitations on the available scan width to avoid the grating lobe formation. This will be taken up a little later.

For attaining the peak value for the antenna main lobe, it is clear that sin(δ) = 0, making the expression using the following L’Hospital’s rule↗.

Hence, sin (δ/2) = sin [(π/λ)×d×(u-u0)] = 0 =sin (±mπ), where m=0,1,2,3,………

So,

or,

In the above equations, parameters u and u0 can have both positive and negative values, but their maximum magnitude is 1 = (|sin(±π/2)|, which is the visible region in sin(θ) space).

Hence,  

Let us now examine as to how this important relationship helps in determining the presence of grating lobes and their angular locations.

  • Case (1): d=λ/2; u= u0±m(λ/d) = u0 ±m×2; u0 is chosen in all these cases to be zero for this study, hence it could be ignored and the equation can be simplified to u= ±m (λ/d), with the condition |±m (λ/d) | max =1.
    • For m=0, u= 0, indicating the broadside beam at θ=0°.
    • For m=1, u=±2 and this not in the visible region.
    • For still higher values of m, the results are not going yield any beam in the visible region. Hence the conclusion is that there is one unique peak and no other lobes of near equal magnitude in the visible range. Hence no grating lobe. This will be the case for all d ≤ λ/2.
  • Case(2): d=λ; u= ±m (λ/d)= ±m;
    • For m=0, u=0; beam at θ=0°, as before.
    • For m=1, u=±1, grating lobes at u= sin (π/2) and u=sin (-π/2)
    • For m =2, u=±2, (not valid for visible range) and this is true for all values of m≥2
    • In conclusion we see, that increase in the inter-distance element d =λ, produces two grating lobes at the edge of the scan range.
  • Case (3): d=1.5 λ; u= ±m (λ/d) = ± 0.66m.
    • For m=0, u=0; beam at θ=0°, as before
    • For m=1, u=± 0.66, grating lobes at θ = sin-1 (± 0.66) = ± 41.3°
    • For m =2, u=±1.34, (not in the visible range; true for all values of m ≥ 2)
  • Case(4): d=2λ; u= ±m (λ/d)= ± 0.5 m;
    • For m=0, u=0, beam at θ=0°, as before.
    • For m=1, u= ± 0.5, so two grating lobes at θ= sin -1 (± 0.5) = ± 30°.
    • For m=2, u= ± 1, indicating another pair of grating lobes at θ=±90°.
    • For m=3, u= ± 1.5 (not valid for visible range). This is the result for all values of m≥ 3. The final tally in this case is the presence of 4 grating lobes along with the main beam
  • Summary
d
(with u0 = 0)
M (Spatial Harmonics)Grating Lobes
≤ λ/2m=0Nil (Main beam only)
λm=0
m=1
Nil(Main beam only);
Present @ u=± 1
1.5λm=0
m=1
Nil(Main beam only);
Present @ u=± 0.66
m=0
m=1
m=2
Nil (Main beam only);
Present  @u=±0.5
Present  @ u=± 1
Table No.1: Inter-Element Spacing Vs. Presence of Grating Lobes

Simulation and results

A simulation run with N=100 in the MATLAB application environment produced the following plot (Fig.2). The plot line in blue indicates the array pattern with d=λ/2 and the line plotted with red color is for the case of d=λ. The third plot in green is for d=1.5 λ. Conclusion arrived in the calculation above are duly verified by this simulation.

Fig.2 Array Pattern for increasing d/λ

The result of a simulation, done to validate the finding in case 4, is given below (Fig.3). The case for d= λ/2 is again taken as reference and blue color line plot indicates its array pattern. The plot with d=2 λ, shows the main beam along with a pair of grating lobes at u=±0.5. Another pair appears u=± 1 (at the edges of the scan) as predicted by the calculation above.

Fig.3: Antenna Pattern for d=λ/2 and 2λ

It will be observed that up to d ≤ λ/2, the higher values of spatial harmonic multiplier (m> 0) is not used, since the value of λ/d 2. It is only the higher values of d/ λ which become smaller fractions on inverting to λ/d, need this harmonic multiplier to satisfy the relationship and answer to the allowable magnitude of |u|max=1. Hence the presence of grating lobes for larger values of d> λ/2.

Designing an array for maximum performance

Having completed this trial, let us now examine the consequence of designing an array for maximum performance at a directed angle of u0=0.5.

The conditional statement at equation (4) is modified to include the requirement of u0 to be set for a value of 0.5:

  • Case 1: Checking the case for λ/d=2 to set the reference: u= -1 ≤ [0.5 ±m (2)] ≤ 1
    • For m=0, u=0.5, giving the required main beam at the set value of 0.5.
    • No other grating lobe with higher harmonic number (m>0) is possible here.
    • So up to d /λ ≤ 0.5, no grating lobes occur as seen earlier.
  • Case 2:For d= λ, u= -1 ≤ [0.5 ±m×(1)] ≤ 1
    • For m=0 : u=0.5 again;
    • For m=1, u=0.5± 1 which yields a grating lobe at u=-0.5; Values higher than m≥2, do not reveal any presence of grating lobe in the visible region.   
  • Case 3: For d=1.5 λ, u= -1 ≤ [0.5 ±m×(0.66)] ≤ 1;
    • For m=0 : u=0.5;
    • For m=1, u= ±0.66, show a pair of grating lobe at θ= ±41.3°.
    • For m=2, u=0.5-1.32=-0.82 yield a grating lobe at θ= -55.1°.
    • Values higher than m≥3, do not reveal any presence of grating lobe in the visible region.
  • Case 4:  For d=2 λ, u= -1 ≤ [0.5 ±m×(0.5)] ≤ 1;
    • For m =0: u=0.5;
    • For m =1, u= {0, 1}, shows a pair of grating lobes at θ=0° and 90°
    • For m=2, u=0.5-1.00=-0.5 yield a grating lobe at θ= -30°.
    • For m=3, u=0.5-1.5 =-1 show a grating lobe at θ= -90°.
    • Once again, values higher than m≥3, do not reveal any presence of grating lobe in the visible region.            
  • Simulation and Results

(i) Occurrence of Grating Lobes:  when d= λ ,(with d= λ/2 as reference)                                          

Fig. 4 Grating lobe for d=λ

The above plot (Fig.4) shows the simulation results for the two cases of d=λ/2 and λ. The validation can be checked through two methods: By plotting the phase data against u, we are able to identify the zeros of the function at which points zero- crossings occur. It will also identify the multiple points where it occurs for d > λ/2. By co-locating the array pattern below, one can easily correlate position of zeros with that of the locations of the grating lobes.

(ii) Similarly, Fig.5 below represents the simulation results and validation for the case d=1.5λ and .

Fig. 5: Grating lobes for d= 1.5λ & 2λ

Limits of Scan imposed by Grating Lobes

Revisiting the relationship for the array to perform in the visible region of u=sin(θ) space, we have

For the safe scanning of the array without the presence of grating lobes, the above condition must be sustained. Let us see how this limit plays out for various values of d (the inter-element spacing), subjected to the limit |u|max=1 for visual range.

It is obvious that as the inter-element space d increases from λ/2, the limits on the permissible scan range gets reduced, if one were to avoid the presence of grating lobes in the operating region. Table 2 compares the allowable scan limits for various values of the inter-element spacing with the constraint that |u|max=1.

d/λElectronic Scanning Limits
(|u0| =0 )
Electronic Scanning Limits
( |u0|>0)
1/2-1 ≤ u ≤ 1 (Full visual range)  -1  ± |u0 | ≤  u ≤  1 ± |u0|
1     -0.5  ≤  u  ≤  0.5    0.5 ± |u0 |  ≤  u ≤  0.5 ± |u0 |
1.5  -0.33  ≤  u  ≤  0.33    -0.33 ± |u0 |  ≤  u ≤  0.33 ± |u0 |  
2  -0.25  ≤  u  ≤  0.25    -0.25 ± |u0 |  ≤  u ≤  0.25  ± |u0 |  
Table 2 : d/λ Vs. Electronic Scan Limits

Simulation results to confirm these limits are presented below for the case of a linear array of 100 elements, with d=1.5, with two conditions: |u0|= 0 and 0.2.

 Case 1: With the condition u0=0, the allowable scan limits are u= ± 0.33;

In the first simulation run, the array scan was arranged to be within allowable limits of u=±0.33 with u0 set to zero. As predicted by the calculation, there is no grating lobe under this situation (Fig. 6). If now, the scan is increased to the full range of visible region (u= ± 1), grating lobes @ u=±o.66, as predicted, appear on either side of the main coverage (Fig. 7) :

Absence of Grating Lobes when scan is within u= ± 0.33; (N=100, d=1.5λ, u0=0)

Fig. 6: No grating lobe for scan range u = ± 0.33

Occurrence of Grating lobes when scan range exceeds limit of u= ± 0.33

Fig. 7 Occurrence of grating lobes when scan range extends to u=± 1

Case 2 : (i) In this simulation, ‘u0‘ is set to a value 0f 0.2 and the simulation repeated as done for the case 1. In the following instance, the permissible range for u is -0.33+0.2 ≤ u ≤ 0.33+0.2; Hence no grating lobe is seen within this region [-0.13≤ u ≤ 0.53].(Fig. 8, below).The top plot gives the locations of the zeros of the phase function. Correspondence of the location of the peak at u=0.2 can be verified with both the plots.

(ii) In the case of the array being scanned to the limits -1≤ u ≤ 1, the appearance of the grating lobes and their angular positions are validated with the calculations shown earlier : (lower plot of Fig.8)

Fig. 8: Grating lobes : When Scan Range extends to u=± 1

Note: To make the scanning range flexible in the simulation, a window function was generated to move along the ‘u’ axis with variable range limits. The following Fig.9, illustrates the results obtained.

A phase function window was generated for a scan range of u= ± 0.33 and ± 0.7. (Blue and Red lines respectively). In the first case the array exhibits no grating lobe as per the calculation. The second case is beyond the permissible scan range and hence shows the presence of grating lobes. Such windowing function helps to evaluate different scan ranges quickly for a designer.

Fig. 9: Grating Lobes (by using window function for ‘u’)

In conclusion, it is stressed that electronic scanning is not without limitation (true for any technology), but it is amenable to system optimization with many of the parameters under the designer’s dynamic control.

Continue reading on Array pattern multiplication

Note: There is a rating embedded within this post, please visit this post to rate it.

Articles in this series

Electronic Scanning Arrays

Articles in this series

Electronic scanning arrays

This is considered to be the best feature of a phased array antenna. To be able to scan in a given space, without having to rotate the antenna is called inertia less scanning and serves many advantages, especially in field operations.

To illustrate the principle, let us take a simple three element linear array arranged with the following phase relationship due to its beam tilted towards an angle θ to receive signal, as shown below:

Figure 1: A simple three element linear array

The path difference between adjacent elements is= d sin (θ) and it progressively increases from right to left. So BB’= d sin (θ) whereas CC’= 2d sin (θ). If we assume that an incoming field with planar wavefront having equal amplitude and equiphase is across CC’, one can visualize the extra path that has to be traveled to reach element A and B, in comparison to the element C. Translating this into phase shift, we have a phase delay= φB = (2π/λ) d sin (θ) at element B, and φA= (2π/λ)(2d)sin(θ) at element A, all reckoned with element C as the point of phase reference. So in the final tally, the field captured along the three-element array comes out to be EResultant = E (1+e-jφ+ e-j2φ), where φ = (2π/λ) d sin(θ) and E is the field value of the incoming equiphase front. The magnitude of the resultant field can be simplified to E(1+2 cos (φ)). It is clear that the field captured by all elements adds up to 3E only when cos (φ) equals zero, which is the case when the antenna is looking at its perpendicular direction, called the broadside. But this is the initial position without scanning. The expression further shows that array gain falls as the angle of scan increases, theoretically falling to E at φ=± π/2. (A simple vector diagram will reveal that at this state, the field received by element A will oppose that of element C). This position also coincides with the limits of scan in the visible region corresponding to a scan limit of θ= ± 90°, for (d/λ) =1/2.

The above example forms the basis for exploring the method of electronic scanning, further.

If we need to receive a signal from an angle inclined from the bore-sight (θ=0°), then there is an inevitable phase shift (delay) to signal reaching the different elements. Hence the uniform phase front along the linear antenna array is no longer along its axis but lifted (refer line CC’) at an angle equal to the angle of scan, (θ in the above example). If we wish to bring back this tilted equiphase plane along the axis of the array, it is obvious each element should now have a phase shift of its own to counter the phase delay they suffered. The required phase shift would be equal in magnitude but in opposite direction.

Let the array be now modified with each element provided with a phase shifter individually. A phase shift + ψ is introduced in a progressive manner from left to right to compensate the phase delays produced earlier by the incoming signal.

Figure 2: Array elements with phase shifters

The modified array equation gives the resultant field = E (1+e-j(φ-ψ) + e-j2(φ-ψ))

Making (φ-ψ) = δ, the above reduces to Ee-j δ(e+jδ + 1 + e-jδ) = E e-j δ [1+2 cos (δ)], which attains maximum at δ =0, making φ=ψ. Remember that φ=(2π/λ) d sin (θ).

Suppose we want the beam to position at θ=30°, then ψ = (2π/λ) d sin (30°). If we need to scan to a direction θ=-30°, then  ψ =- (2π/λ) d sin(30°), will be a lagging phase shift.

What essentially resulted was that the equiphase plane inclined at θ was brought back to the array surface by phase weighting (ψ) of the element as arranged in Figure 2. The phase weight is dependent on the angle chosen for the scan, once the designer has already selected the spacing among elements (d/λ). This phase weight will then have to continually track its value in sync with the position of the scan angle ‘θ ‘.

A general equation can now be attempted to describe the pattern of a linear array antenna with N elements:

A simulation run with an array antenna is shown in Figure 1. A linear array of 100 elements was chosen (N=100), with the inter-distance between elements d=λ/2. The desired angle of scan was u1=sin(θ0) This was set for three values: u=0 (broadside), and at u1=± 0.5 (angle of ±30°).

Figure 3: Array scan pattern

It is true that the gain of the scanned patterns becomes a function of the scanned angle [Gθ = G cos(θ)]; this fact is not visible in this plot as each pattern is normalized with respect to its own peak.

Continue reading about the discussion on grating lobes…

For further reading

Benjamin et. al., “Phased array technologies”, Proceedings of the IEEE, Vol. 104, No. 3, March 201

Rate this post: Note: There is a rating embedded within this post, please visit this post to rate it.

Phased Array Antenna – an introduction

Key focus: Phased Array Antenna, employs simpler elemental antennas to configure a complex system, has many interesting features and capabilities. 

Articles in this series

Since the latter half of previous century, a class of antenna technology called the ‘Phased Arrays’ has witnessed phenomenal progress in design, engineering and applications. Concurrent progress in solid state material has enabled active antenna elements to be packed as arrays in many different geometries- linear, planar and conformal shapes. Though started as the solution for various taxing military needs, this technology has since entered many civilian domains, the space applications being one of the prominent usage in the present times.

Phased Array, which employs simpler elemental antennas to configure a complex system, has many interesting features and capabilities. The first one that appeals to a designer, is the ability to address the design problem by reaching at the elemental level for control. Theoretically, it is possible to control at each element, the amplitude and phase of the signal being fed. This gives rise to a dynamic control of antenna beam shape, electronic scanning, sidelobes and their placement-all leading to a variety of applications from the same system.

Though amplitude and phase control are possible at element level, it is generally seen that the designers use the amplitude control to set the illumination function of the elements and deal with the dynamic requirements by phase shifter switching. Presently, phase shifters are capable of being engineered in small form factor, incorporating the solid state source, the connecting geometry and the radiating element. Digital and modular technology make it feasible for mass production and for maintainability and reliability in severe environmental condition.

The phase shifters employed in such systems use solid state devices and electronically controlled ferrites for switching.  They are typically controlled with 4 to 5 bits digital accuracy which give enough phase states to design for. The setting accuracy of a phase state has also considerably advanced in that a designer can plan for dynamic beam shaping with confidence. By dividing the entire R.F. source among many elements, the power management and heat dissipation of the system get distributed to manageable engineering practice. Hence such a level of sophistication on phase control has finally evolved an array system, justifiably called the ‘Phased Array Antenna’.

In the succeeding pages, effort has been made to describe what is the unique specialty of Phased Arrays- ‘the electronic steering’ –its value and also the limitations, so that a balanced view is derived. The idea is to employ simple examples to start with and build on this initial strength to move onto more intricate antenna functions.  The emphasis will be on the design concepts and their practical utility. No deeper background in EM theory, RF-Microwave practices are required at present; but essential knowledge of antennas and their different operating frequency environments are already visible to everyone in today’s technology bound world.

It is presumed that readers are familiar with many of the antenna concepts and parameters, like antenna gain, beamwidth, sidelobes, field and power patterns, reciprocity in antenna for transmitting and receiving; and allied engineering mathematical concepts like Fourier Transform, geometric series and their summation, and the need for logarithmic compression in the plotting of antenna patterns. Such an initial preparation will help in assimilating the text and graphics which follow.

Continue reading the discussion on electronic scanning arrays… (not yet published! Monitor this place for the follow-up article)

Rate this post: Note: There is a rating embedded within this post, please visit this post to rate it.

Chirp Signal – FFT & PSD in Matlab & Python

Key focus: Know how to generate a Chirp signal, compute its Fourier Transform using FFT and power spectral density (PSD) in Matlab & Python.

This article is part of the following books
Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885
Digital Modulations using Python ISBN: 978-1712321638
Wireless communication systems in Matlab ISBN: 979-8648350779
All books available in ebook (PDF) and Paperback formats

Introduction

All the signals discussed so far do not change in frequency over time. Obtaining a signal with time-varying frequency is of main focus here. A signal that varies in frequency over time is called “chirp”. The frequency of the chirp signal can vary from low to high frequency (up-chirp) or from high to low frequency (low-chirp).

Observation

Chirp signals/signatures are encountered in many applications ranging from radar, sonar, spread spectrum, optical communication, image processing, doppler effect, motion of a pendulum, as gravitation waves, manifestation as Frequency Modulation (FM), echo location [1] etc.

Mathematical Description:

A linear chirp signal sweeps the frequency from low to high frequency (or vice-versa) linearly. One approach to generate a chirp signal is to concatenate a series of segments of sine waves each with increasing(or decreasing) frequency in order. This method introduces discontinuities in the chirp signal due to the mismatch in the phases of each such segments. Modifying the equation of a sinusoid to generate a chirp signal is a better approach.

The equation for generating a sinusoidal (cosine here) signal with amplitude A, angular frequency and initial phase is

This can be written as a function of instantaneous phase

where is the instantaneous phase of the sinusoid and it is linear in time. The time derivative of instantaneous phase is equal to the angular frequency of the sinusoid – which in case is a constant in the above equation.

Instead of having the phase linear in time, let’s change the phase to quadratic form and thus non-linear.

for some constant .

Therefore, the equation for chirp signal takes the following form,

The first derivative of the phase, which is the instantaneous angular frequency becomes a function of time, which is given by

The time-varying frequency in Hertz is given by

In the above equation, the frequency is no longer a constant, rather it is of time-varying nature with initial frequency given by . Thus, from the above equation, given a time duration T, the rate of change of frequency is given by

where, is the starting frequency of the sweep, is the frequency at the end of the duration T.

Substituting (7) & (8) in (6)

From (6) and (8)

where  is a constant which will act as the initial phase of the sweep.

Thus the modified equation for generating a chirp signal (from equations (5) and (10)) is given by

where the time-varying frequency function is given by

Generation of Chirp signal, computing its Fourier Transform using FFT and power spectral density (PSD) in Matlab is shown as example, for Python code, please refer the book Digital Modulations using Python.

Generating a chirp signal without using in-built “chirp” Function in Matlab:

Implement a function that describes the chirp using equation (11) and (12). The starting frequency of the sweep is and the frequency at time is . The initial phase forms the final part of the argument in the following function

function x=mychirp(t,f0,t1,f1,phase)
%Y = mychirp(t,f0,t1,f1) generates samples of a linear swept-frequency
%   signal at the time instances defined in timebase array t.  The instantaneous
%   frequency at time 0 is f0 Hertz.  The instantaneous frequency f1
%   is achieved at time t1.
%   The argument 'phase' is optional. It defines the initial phase of the
%   signal degined in radians. By default phase=0 radian
    
if nargin==4
    phase=0;
end
    t0=t(1);
    T=t1-t0;
    k=(f1-f0)/T;
    x=cos(2*pi*(k/2*t+f0).*t+phase);
end

The following wrapper script utilizes the above function and generates a chirp with starting frequency at the start of the time base and at which is the end of the time base. From the PSD plot, it can be ascertained that the signal energy is concentrated only upto 25 Hz

fs=500; %sampling frequency
t=0:1/fs:1; %time base - upto 1 second

f0=1;% starting frequency of the chirp
f1=fs/20; %frequency of the chirp at t1=1 second
x = mychirp(t,f0,1,f1); 
subplot(2,2,1)
plot(t,x,'k');
title(['Chirp Signal']);
xlabel('Time(s)');
ylabel('Amplitude');

FFT and power spectral density

As with other signals, describes in the previous posts, let’s plot the FFT of the generated chirp signal and its power spectral density (PSD).

L=length(x);
NFFT = 1024;
X = fftshift(fft(x,NFFT));
Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
f = fs*(-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector

subplot(2,2,2)
plot(f,abs(X)/(L),'r');
title('Magnitude of FFT');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');
xlim([-50 50])


Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
subplot(2,2,3)
plot(f,10*log10(Pxx),'r');
title('Double Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
xlim([-100 100])

X = fft(x,NFFT);
X = X(1:NFFT/2+1);%Throw the samples after NFFT/2 for single sided plot
Pxx=X.*conj(X)/(NFFT*NFFT);
f = fs*(0:NFFT/2)/NFFT; %Frequency Vector
subplot(2,2,4)
plot(f,10*log10(Pxx),'r');
title('Single Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');

For Python code, please refer the book Digital Modulations using Python

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

References:

[1] Patrick Flandrin,“Chirps everywhere”,CNRS — Ecole Normale Supérieure de Lyon

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Gaussian Pulse – FFT & PSD in Matlab & Python

Key focus: Know how to generate a gaussian pulse, compute its Fourier Transform using FFT and power spectral density (PSD) in Matlab & Python.

Numerous texts are available to explain the basics of Discrete Fourier Transform and its very efficient implementation – Fast Fourier Transform (FFT).  Often we are confronted with the need to generate simple, standard signals (sine, cosine, Gaussian pulse, squarewave, isolated rectangular pulse, exponential decay, chirp signal) for simulation purpose. I intend to show (in a series of articles) how these basic signals can be generated in Matlab and how to represent them in frequency domain using FFT.

This article is part of the following books
Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885
Digital Modulations using Python ISBN: 978-1712321638
Wireless communication systems in Matlab ISBN: 979-8648350779
All books available in ebook (PDF) and Paperback formats

Gaussian Pulse : Mathematical description:

In digital communications, Gaussian Filters are employed in Gaussian Minimum Shift Keying – GMSK (used in GSM technology) and Gaussian Frequency Shift Keying (GFSK). Two dimensional Gaussian Filters are used in Image processing to produce Gaussian blurs. The impulse response of a Gaussian Filter is Gaussian. Gaussian Filters give no overshoot with minimal rise and fall time when excited with a step function. Gaussian Filter has minimum group delay. The impulse response of a Gaussian Filter is written as a Gaussian Function as follows

The Fourier Transform of a Gaussian pulse preserves its shape.

The above derivation makes use of the following result from complex analysis theory and the property of Gaussian function – total area under Gaussian function integrates to 1.

By change of variable, let ( ). 

Thus, the Fourier Transform of a Gaussian pulse is a Gaussian Pulse.

Gaussian Pulse – Fourier Transform using FFT (Matlab & Python):

The following code generates a Gaussian Pulse with ( ). The Discrete Fourier Transform of this digitized version of Gaussian Pulse is plotted with the help of (FFT) function in Matlab.

For Python code, please refer the book Digital Modulations using Python

fs=80; %sampling frequency
sigma=0.1;
t=-0.5:1/fs:0.5; %time base

variance=sigma^2;
x=1/(sqrt(2*pi*variance))*(exp(-t.^2/(2*variance)));
subplot(2,1,1)
plot(t,x,'b');
title(['Gaussian Pulse \sigma=', num2str(sigma),'s']);
xlabel('Time(s)');
ylabel('Amplitude');

L=length(x);
NFFT = 1024;
X = fftshift(fft(x,NFFT));
Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
f = fs*(-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector

subplot(2,1,2)
plot(f,abs(X)/fs,'r');
title('Magnitude of FFT');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');
xlim([-10 10])
Figure 1: Gaussian pulse and its FFT (magnitude)

Double Sided and Single Power Spectral Density using FFT:

Next, the Power Spectral Density (PSD) of the Gaussian pulse is constructed using the FFT. PSD describes the power contained at each frequency component of the given signal. Double Sided power spectral density is plotted first, followed by single sided power spectral density plot (retaining only the positive frequency side of the spectrum).

Pxx=X.*conj(X)/(L*L); %computing power with proper scaling
figure;
plot(f,10*log10(Pxx),'r');
title('Double Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
Figure 2: Double sided power spectral density of Gaussian pulse
X = fft(x,NFFT);
X = X(1:NFFT/2+1);%Throw the samples after NFFT/2 for single sided plot
Pxx=X.*conj(X)/(L*L);
f = fs*(0:NFFT/2)/NFFT; %Frequency Vector
figure;
plot(f,10*log10(Pxx),'r');
title('Single Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');

For Python code, please refer the book Digital Modulations using Python

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Generating Basic signals – Rectangular Pulse and Power Spectral Density using FFT

Numerous texts are available to explain the basics of Discrete Fourier Transform and its very efficient implementation – Fast Fourier Transform (FFT).  Often we are confronted with the need to generate simple, standard signals (sine, cosineGaussian pulsesquare waveisolated rectangular pulse, exponential decay, chirp signal) for simulation purpose. I intend to show (in a series of articles) how these basic signals can be generated in Matlab and how to represent them in frequency domain using FFT.

This article is part of the book Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885 available in ebook (PDF) format (click here) and Paperback (hardcopy) format (click here)
Wireless Communication Systems in Matlab, ISBN: 978-1720114352 available in ebook (PDF) format (click here) and Paperback (hardcopy) format (click here).

Rectangular pulse: mathematical description

An isolated rectangular pulse of amplitude A and duration T is represented mathematically as

where

The Fourier transform of isolated rectangular pulse g(t) is

where, the sinc function is given by

Thus, the Fourier Transform pairs are

The Fourier Transform describes the spectral content of the signal at various frequencies. For a given signal g(t), the Fourier Transform is given by

where, the absolute value gives the magnitude of the frequency components (amplitude spectrum) and are their corresponding phase (phase spectrum) . For the rectangular pulse, the amplitude spectrum is given as

The amplitude spectrum peaks at f=0 with value equal to AT. The nulls of the spectrum occur at integral multiples of 1/T, i.e, ( )

Generating an isolated rectangular pulse in Matlab:

An isolated rectangular pulse of unit amplitude and width w (the factor T in equations above ) can be generated easily with the help of in-built function – rectpuls(t,w) command in Matlab. As an example, a unit amplitude rectangular pulse of duration is generated.

fs=500; %sampling frequency
T=0.2; %width of the rectangule pulse in seconds

t=-0.5:1/fs:0.5; %time base

x=rectpuls(t,T); %generating the square wave

plot(t,x,'k');
title(['Rectangular Pulse width=', num2str(T),'s']);
xlabel('Time(s)');
ylabel('Amplitude');

Amplitude spectrum using FFT:

Matlab’s FFT function is utilized for computing the Discrete Fourier Transform (DFT). The magnitude of FFT is plotted. From the following plot, it can be noted that the amplitude of the peak occurs at f=0 with peak value  . The nulls in the spectrum are located at  ().

L=length(x);
NFFT = 1024;
X = fftshift(fft(x,NFFT)); %FFT with FFTshift for both negative & positive frequencies
f = fs*(-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector

figure;
plot(f,abs(X)/(L),'r');
title('Magnitude of FFT');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');

Power spectral density (PSD) using FFT:

The distribution of power among various frequency components is plotted next. The first plot shows the double-side Power Spectral Density which includes both positive and negative frequency axis. The second plot describes the PSD only for positive frequency axis (as the response is just the mirror image of negative frequency axis).

figure;
Pxx=X.*conj(X)/(L*L); %computing power with proper scaling
plot(f,10*log10(Pxx),'r');
title('Double Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
X = fft(x,NFFT);
X = X(1:NFFT/2+1);%Throw the samples after NFFT/2 for single sided plot
Pxx=X.*conj(X)/(L*L);
f = fs*(0:NFFT/2)/NFFT; %Frequency Vector
plot(f,10*log10(Pxx),'r');
title('Single Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');

Magnitude and phase spectrum:

The phase spectrum of the rectangular pulse manifests as series of pulse trains bounded between 0 and , provided the rectangular pulse is symmetrically centered around sample zero. This is explained in the reference here and the demo below.

clearvars;
x = [ones(1,7) zeros(1,127-13) ones(1,6)];
subplot(3,1,1); plot(x,'k');
title('Rectangular Pulse'); xlabel('Sample#'); ylabel('Amplitude');

NFFT = 127;
X = fftshift(fft(x,NFFT)); %FFT with FFTshift for both negative & positive frequencies
f = (-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector

subplot(3,1,2); plot(f,abs(X),'r');
title('Magnitude Spectrum'); xlabel('Frequency (Hz)'); ylabel('|X(f)|');

subplot(3,1,3); plot(f,atan2(imag(X),real(X)),'r');
title('Phase Spectrum'); xlabel('Frequency (Hz)'); ylabel('\angle X(f)');

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Generating Basic signals – Square Wave and Power Spectral Density using FFT

Note: There is a rating embedded within this post, please visit this post to rate it.

Numerous texts are available to explain the basics of Discrete Fourier Transform and its very efficient implementation – Fast Fourier Transform (FFT).  Often we are confronted with the need to generate simple, standard signals (sine, cosineGaussian pulsesquarewaveisolated rectangular pulse, exponential decay, chirp signal) for simulation purpose. I intend to show (in a series of articles) how these basic signals can be generated in Matlab and how to represent them in frequency domain using FFT.

This article is part of the book Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885 available in ebook (PDF) format (click here) and Paperback (hardcopy) format (click here)
Wireless Communication Systems in Matlab, ISBN: 978-1720114352 available in ebook (PDF) format (click here) and Paperback (hardcopy) format (click here).

Significance of Square Waves

The most logical way of transmitting information across a communication channel is through a stream of square pulse – a distinct pulse for ‘0‘ and another for ‘1‘. Digital signals are graphically represented as square waves with certain symbol/bit period. Square waves are also used universally in switching circuits, as clock signals synchronizing various blocks of digital circuits, as reference clock for a given system domain and so on.

Square wave manifests itself as a wide range of harmonics in frequency domain and therefore can cause electromagnetic interference. Square waves are periodic and contain odd harmonics when expanded as Fourier Series (where as signals like saw-tooth and other real word signals contain harmonics at all integer frequencies). Since a square wave literally expands to infinite number of odd harmonic terms in frequency domain, approximation of square wave is another area of interest. The number of terms of its Fourier Series expansion, taken for approximating the square wave is often seen as Gibbs Phenomenon, which manifests as ringing effect at the corners of the square wave in time domain (visual explanation here).

True Square waves are a special class of rectangular waves with 50% duty cycle. Varying the duty cycle of a rectangular wave leads to pulse width modulation, where the information is conveyed by changing the duty-cycle of each transmitted rectangular wave.

How to generate a square wave in Matlab

If you know the trick of generating a sine wave in Matlab, the task is pretty much simple. Square wave is generated using “square” function in Matlab. The command sytax – square(t,dutyCycle) – generates a square wave with period for the given time base. The command behaves similar to “sin” command (used for generating sine waves), but in this case it generates a square wave instead of a sine wave. The argument – dutyCycle is optional and it defines the desired duty cycle of the square wave. By default (when the dutyCycle argument is not supplied) the square wave is generated with (50%) duty cycle.

f=10; %frequency of sine wave in Hz
overSampRate=30; %oversampling rate
fs=overSampRate*f; %sampling frequency
duty_cycle=50; % Square wave with 50% Duty cycle (default)
nCyl = 5; %to generate five cycles of sine wave

t=0:1/fs:nCyl*1/f; %time base

x=square(2*pi*f*t,duty_cycle); %generating the square wave

plot(t,x,'k');
title(['Square Wave f=', num2str(f), 'Hz']);
xlabel('Time(s)');
ylabel('Amplitude');

Power Spectral Density using FFT

Let’s check out how the generated square wave will look in frequency domain. The Fast Fourier Transform (FFT) is utilized here. As discussed in the article here, there are numerous ways to plot the response of FFT. Single Sided power spectral density is plotted first, followed by the Double-sided power spectral density.

Single Sided Power Spectral Density

X = fft(x,NFFT);
X = X(1:NFFT/2+1);%Throw the samples after NFFT/2 for single sided plot
Pxx=X.*conj(X)/(NFFT*L);
f = fs*(0:NFFT/2)/NFFT; %Frequency Vector

plot(f,10*log10(Pxx),'r');
title('Single Sided Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
ylim([-45 -5])

Double Sided Power Spectral Density

L=length(x);
NFFT = 1024;
X = fftshift(fft(x,NFFT));
Pxx=X.*conj(X)/(NFFT*L); %computing power with proper scaling
f = fs*(-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector

plot(f,10*log10(Pxx),'r');
title('Double Sided Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
Note: There is a rating embedded within this post, please visit this post to rate it.

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Plot FFT using Matlab – FFT of sine wave & cosine wave

Key focus: Learn how to plot FFT of sine wave and cosine wave using Matlab. Understand FFTshift. Plot one-sided, double-sided and normalized spectrum.

Introduction

Numerous texts are available to explain the basics of Discrete Fourier Transform and its very efficient implementation – Fast Fourier Transform (FFT).  Often we are confronted with the need to generate simple, standard signals (sine, cosine, Gaussian pulse, squarewave, isolated rectangular pulse, exponential decay, chirp signal) for simulation purpose. I intend to show (in a series of articles) how these basic signals can be generated in Matlab and how to represent them in frequency domain using FFT. If you are inclined towards Python programming, visit here.

This article is part of the following books
Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885
Wireless communication systems in Matlab ISBN: 979-8648350779
All books available in ebook (PDF) and Paperback formats

Sine Wave

In order to generate a sine wave in Matlab, the first step is to fix the frequency of the sine wave. For example, I intend to generate a f=10 Hz sine wave whose minimum and maximum amplitudes are and respectively. Now that you have determined the frequency of the sinewave, the next step is to determine the sampling rate. Matlab is a software that processes everything in digital. In order to generate/plot a smooth sine wave, the sampling rate must be far higher than the prescribed minimum required sampling rate which is at least twice the frequency – as per Nyquist Shannon Theorem. A oversampling factor of is chosen here – this is to plot a smooth continuous-like sine wave (If this is not the requirement, reduce the oversampling factor to desired level). Thus the sampling rate becomes . If a phase shift is desired for the sine wave, specify it too.

f=10; %frequency of sine wave
overSampRate=30; %oversampling rate
fs=overSampRate*f; %sampling frequency
phase = 1/3*pi; %desired phase shift in radians
nCyl = 5; %to generate five cycles of sine wave

t=0:1/fs:nCyl*1/f; %time base

x=sin(2*pi*f*t+phase); %replace with cos if a cosine wave is desired
plot(t,x);
title(['Sine Wave f=', num2str(f), 'Hz']);
xlabel('Time(s)');
ylabel('Amplitude');

Representing in Frequency Domain

Representing the given signal in frequency domain is done via Fast Fourier Transform (FFT) which implements Discrete Fourier Transform (DFT) in an efficient manner. Usually, power spectrum is desired for analysis in frequency domain. In a power spectrum, power of each frequency component of the given signal is plotted against their respective frequency. The command computes the -point DFT. The number of points – –  in the DFT computation is taken as power of (2) for facilitating efficient computation with FFT. A value of is chosen here. It can also be chosen as next power of 2 of the length of the signal.

Different representations of FFT:

Since FFT is just a numeric computation of -point DFT, there are many ways to plot the result.

1. Plotting raw values of DFT:

The x-axis runs from to – representing sample values. Since the DFT values are complex, the magnitude of the DFT is plotted on the y-axis. From this plot we cannot identify the frequency of the sinusoid that was generated.

NFFT=1024; %NFFT-point DFT      
X=fft(x,NFFT); %compute DFT using FFT        
nVals=0:NFFT-1; %DFT Sample points       
plot(nVals,abs(X));      
title('Double Sided FFT - without FFTShift');        
xlabel('Sample points (N-point DFT)')        
ylabel('DFT Values');

2. FFT plot – plotting raw values against Normalized Frequency axis:

In the next version of plot, the frequency axis (x-axis) is normalized to unity. Just divide the sample index on the x-axis by the length of the FFT. This normalizes the x-axis with respect to the sampling rate . Still, we cannot figure out the frequency of the sinusoid from the plot.

NFFT=1024; %NFFT-point DFT      
X=fft(x,NFFT); %compute DFT using FFT        
nVals=(0:NFFT-1)/NFFT; %Normalized DFT Sample points         
plot(nVals,abs(X));      
title('Double Sided FFT - without FFTShift');        
xlabel('Normalized Frequency')       
ylabel('DFT Values');

3. FFT plot – plotting raw values against normalized frequency (positive & negative frequencies):

As you know, in the frequency domain, the values take up both positive and negative frequency axis. In order to plot the DFT values on a frequency axis with both positive and negative values, the DFT value at sample index has to be centered at the middle of the array. This is done by using function in Matlab. The x-axis runs from to where the end points are the normalized ‘folding frequencies’ with respect to the sampling rate .

NFFT=1024; %NFFT-point DFT      
X=fftshift(fft(x,NFFT)); %compute DFT using FFT      
fVals=(-NFFT/2:NFFT/2-1)/NFFT; %DFT Sample points        
plot(fVals,abs(X));      
title('Double Sided FFT - with FFTShift');       
xlabel('Normalized Frequency')       
ylabel('DFT Values');

4. FFT plot – Absolute frequency on the x-axis Vs Magnitude on Y-axis:

Here, the normalized frequency axis is just multiplied by the sampling rate. From the plot below we can ascertain that the absolute value of FFT peaks at and . Thus the frequency of the generated sinusoid is . The small side-lobes next to the peak values at and are due to spectral leakage.

NFFT=1024;      
X=fftshift(fft(x,NFFT));         
fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT;        
plot(fVals,abs(X),'b');      
title('Double Sided FFT - with FFTShift');       
xlabel('Frequency (Hz)')         
ylabel('|DFT Values|');

5. Power Spectrum – Absolute frequency on the x-axis Vs Power on Y-axis:

The following is the most important representation of FFT. It plots the power of each frequency component on the y-axis and the frequency on the x-axis. The power can be plotted in linear scale or in log scale. The power of each frequency component is calculated as

Where is the frequency domain representation of the signal . In Matlab, the power has to be calculated with proper scaling terms (since the length of the signal and transform length of FFT may differ from case to case).

NFFT=1024;
L=length(x);         
X=fftshift(fft(x,NFFT));         
Px=X.*conj(X)/(NFFT*L); %Power of each freq components       
fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT;        
plot(fVals,Px,'b');      
title('Power Spectral Density');         
xlabel('Frequency (Hz)')         
ylabel('Power');

If you wish to verify the total power of the signal from time domain and frequency domain plots, follow this link.
Plotting the power spectral density (PSD) plot with y-axis on log scale, produces the most encountered type of PSD plot in signal processing.

NFFT=1024;      
L=length(x);         
X=fftshift(fft(x,NFFT));         
Px=X.*conj(X)/(NFFT*L); %Power of each freq components       
fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT;        
plot(fVals,10*log10(Px),'b');        
title('Power Spectral Density');         
xlabel('Frequency (Hz)')         
ylabel('Power');

6. Power Spectrum – One-Sided frequencies

In this type of plot, the negative frequency part of x-axis is omitted. Only the FFT values corresponding to to sample points of -point DFT are plotted. Correspondingly, the normalized frequency axis runs between to . The absolute frequency (x-axis) runs from to .

L=length(x);        
NFFT=1024;       
X=fft(x,NFFT);       
Px=X.*conj(X)/(NFFT*L); %Power of each freq components       
fVals=fs*(0:NFFT/2-1)/NFFT;      
plot(fVals,Px(1:NFFT/2),'b','LineSmoothing','on','LineWidth',1);         
title('One Sided Power Spectral Density');       
xlabel('Frequency (Hz)')         
ylabel('PSD');

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

For further reading

[1] Power spectral density – MIT opencourse ware↗

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author


Wireless Communication Systems in Matlab
Second Edition(PDF)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Python
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart

Digital Modulations using Matlab
(PDF ebook)

Note: There is a rating embedded within this post, please visit this post to rate it.
Checkout Added to cart
Hand-picked Best books on Communication Engineering
Best books on Signal Processing