From ControlTheoryPro.com

Contents
1 Introduction to psdData Class Example
In Stochastic Controls Power Spectral Densities (PSDs) are often used to represent or define disturbance and noise environments. However, there are several implementation issues that must be dealt with when using PSDs.
The basics of a PSD are as follows:
 The PSD is represented in MATLAB by a vector of frequencies and magnitudes at those frequencies.
 PSDs typically (but not always) have units of magnitude squared (^2) over Hz.
 The primary metric used to compare one PSD to another is the Root Mean Square (RMS).
 As with most metrics the RMS does not tell you everything.
 Particularly important is where in the frequency range does most of the PSD energy resides. More specifically, "Is the energy in the low, mid, or high frequency range by comparison to your designed rejection and performance requirements?"
Some of the PSD implementation issues are:
 PSD measurements are often very coarse particularly in the low frequencies. The result of this is that low frequency measurements might be at 0.01 Hz, 0.1 Hz, 1 Hz, 2 Hz, 5 Hz, 10 Hz, ... For noise PSDs this is not usually an issue. However, disturbance PSDs are often dominated by low frequency content. The RMS is an integration under the curve and coarse measurements often lead to inaccurate RMS calculations.
 Addition and subtraction of PSDs must occur at the same frequency and this often dictates the use of an interpolated PSD.
 Measured data comes in the form that is easiest to capture with the least noise. This often leads to measurements in rate (from good MEMS gyros) or acceleration (from good MEMS accelerometers). While the measurements are in rate or acceleration the disturbances that need to be in any simulation may need to be integrated to a position PSD.
2 Example from PSD Data
In controls for pointing and tracking, the customer's performance requirements often state that the system must "perform such that the residual error (RMS) is no greater than X assuming a disturbance spectrum as seen in Figure Y and Table Z." Let's assume that Table Z is a table of frequencies (in Hz) and magnitudes (in rad) representing an angular position disturbance envelope.
For the purposes of our example let's create some fake angular position PSD data.
Frequency (Hz)  Magnitude (rad^2 / Hz)  

0 
1 

1 
1 

10 
1E3 

100 
1E9 

1000 
1E9 
The resulting MATLAB Code is
>> freq = [0 ; 1 ; 10 ; 100 ; 1000]; >> mag = [1 ; 1 ; 1E3 ; 1E9 ; 1E9];
Example Instantiation:
>> a = psdData('Strawman PSD Data', freq, mag, 'psd', [], 'interpolate', 5000);
Since the data is in rad and the default for the psdData Class is rad/sec we need to change the units
>> setUnits(a, 'Magnitude', 'rad', 'TimeUnit', 'Position')
The command above set the units for the magnitude to rad (radians).
2.1 Plot
Plot Command:
>> plot([], a, 'rms legend', 'southwest')
For more detailed information on the plotting functions see the plot method.
2.2 Differentiate
Differentiate Command:
>> b = differentiate(a);
Plot Command:
>> plot([], b, 'rms legend', 'southwest')
Note that the Name property of the psdData class changes when integrate or differentiate are used to reflect the operation. The units are automatically updated as well. A more meaningful name might be
>> b.Name = 'Strawman Rate PSD Data';
3 Example from Time Series Data
Example Instantiation:
>> a = psdData('Example Time Series Data', time, data, 'time', window, 'interpolate', 5000);
where window determines the size of the window used in creating the PSD and can be any number of seconds up the maximum time of the vector time.
However, note that that the original data had a transient at the very begining. This is likely an artifact of the measurement process and a it would be more representative to create the PSD from the portion of the Time Series data that came after this transient had passed. It appears that the transient is over by second 2 so create a new psdData object like this
>> ind = find(time >= 2); >> newTime = time(ind)  min(time(ind)); >> b = psdData('Trimmed Time Series Data', ... newTime, ... data(ind), ... 'time', ... window, ... 'interpolate', 5000);
The variable ind captures all indices in the time vector that are greater than or equal to 2. However, for the PSD to be calculated correctly the time vector needs to start at 0. So the newTime variable was created to hold a shifted time vector. Now the psdData object will construct a PSD from the trimmed Time Series data set which is uncorrupted by the startup transient.
3.1 Integrate
Let's assume that this measured data is rate data in rad/sec  the psdData Class default. For our imaginary application we need a disturbance in angular position instead of rate. Therefore we need to integrate the PSD so that it represents a position disturbance. With the psdData class the code is simple
>> positionPSD = integrate(b); >> positionPSD.Name = 'Angular Position Disturbance PSD';
The integrate method will automatically generate a name for the integrated PSD which, in this example, would be of the form Integrated Example Time Series Data. However, I chose to rename the new PSD to give it a more meaningful name. This is important mainly in that the psdData Class name is used when the data is plotted.
Also note that the resulting psdData Class object of the integrate method (as well as the overloaded plus & minus and differentiate methods) will no longer have any time series data. There are several reasons for this but the main one is that direct manipulation of test data without any consideration of its source is problematic at best.
4 Download psdData Class
5 Notes
 psdData assumes that any PSD provided to it is 1sided
 psdData creates 1sided PSDs from Time Series Data
 The integrate, differentiate, plus, and minus methods will provide results based on the PSD meaning that Time Series data is not present in the resulting psdData object (it will still be present in the original objects though)
 Calculation of the PSD from time series data requires the pwelch function and pwelch is part of the Signal Processing Toolbox.