How can I generate power spectral density from waveform?

Hello. I would like to generate a waveform in MATLAB according to an equation, and then calculate its power spectral density (PSD). I have some experience with MATLAB, but I admit that my DSP skills need improving, and I think it's the latter which is causing me some troubles...

Waveform generation

The waveform is data that has been modulated by a specific method onto a sinusoidal carrier. The waveform can be represented mathematically by:

where:

• is the RF carrier frequency
• m is the modualtion index, which can be set to 1.2 rad
• is a binary data sequence with symbol rate

depends on the modulation scheme; I give the detail so you can follow the code below:

• For PCM/PM:
• For PCM/PSK/PM with square-wave subcarrier:
• For PCM/PSK/PM with sine-wave subcarrier:

where:

• is the subcarrier frequency
• is the subcarrier phase

I have some code that successfully generates the waveforms:

%% Set up PRBS bit stream, "b". Make function of time after defining samples per bit
clear; close all; clc;
n = 2^6; % # of bits
rng default; % set rng seed to default to create repeatable results
b = randi([0,1],n,1); % set up PRBS between 0 and +1
b(b==0) = -1; % change all 0s to -1 so data is -1 and +1
b = b';
spb = 2^6; % samples per bit
ns = n*spb; % total # of samples
b = repelem(b,spb); % recreate bit stream, but repeating the elements so you have "spb" samples per bit
Rb = 2000; % bit rate = 2kb/s
Tb = 1/Rb; % bit period
Ts = Tb/spb; % sample period
t = 0:Ts:(ns-1)*Ts; % set up time array
plot_lim = spb*6; % used to plot arbitrary portion of bit stream / waveform. i.e. plot 6 bits worth of data
fig = figure; 
subplot(1,2,1); plot(t(1:plot_lim),b(1:plot_lim)); title('Bit stream'); xlabel('t [s]');

% Generate waveforms
f_c = 0; % carrier frequency
omega_c = 2*pi*f_c; 
m = 1.2; % modultaion index
sc_Rb_ratio = 4; % sub-carrier to symbol rate ratio
f_sc = f_c + sc_Rb_ratio*Rb; % sub-carrier frequency
omega_sc = 2*pi*f_sc;
theta_sc = 0; % sub-carrier phase
Pt = [repelem(1,length(t)); square(omega_sc*t + theta_sc); sin(omega_sc*t + theta_sc)];
s = sqrt(2)*sin(omega_c*t + m*Pt.*b);
subplot(1,2,2); 
plot(t(1:plot_lim),s(:,1:plot_lim)); xlabel('t [s]');
title('Waveform');
legArr = {'PCM/PM';'PCM/PSK/PM (square-wave)';'PCM/PSK/PM (sine-wave)'};
legend(legArr);
fig.Position = [200 300 1200 500];

This results in the following bit stream / waveforms:

Power spectral density

How do I work out the PSD (in the frequency domain) from the defined waveforms? I want to be able to select m and , generate the waveform (in the time domain) and see how it changes the PSD (in the frequency domain).

I would appreciate seeing other people's methodologies to get an idea of the best way to do this. FYI I have access to the Commucations Toolbox and Signal Processing Toolbox. Thank you in advance.

[dpb]

Well, you generate the waveform simply by coding the equation in MATLAB syntax -- which is very close to the physical form excepting you need explicit multiplication operators. Set the t vector as desired and use "dot" operators to vectorize, of course.

There's example of how all this is done, plus the PSD in the doc for FFT; there are alternative ways to compute a PSD in both toolboxes; I don't have Comm TB so no knowledge of it; there may be techniques in it more directly applicable to your subject area than the base methods; they'll all be based on same/similar ideas/techniques, though.

[rkoz575]

Thank you for your comment. I have added some code for generating the waveforms, but am still stuck on the PSD question..

[dpb]

What have you tried and what, specifically stumped you?

Answer 1

[rkoz575]

Hi I have found the answer. One should follow the steps laid out in:

https://uk.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html

This yields the expected spectrum.