|
1 | 1 | """
|
2 | 2 | ========
|
3 |
| -Psd Demo |
| 3 | +PSD Demo |
4 | 4 | ========
|
5 | 5 |
|
6 | 6 | Plotting Power Spectral Density (PSD) in Matplotlib.
|
|
9 | 9 | many useful libraries for computing a PSD. Below we demo a few examples
|
10 | 10 | of how this can be accomplished and visualized with Matplotlib.
|
11 | 11 | """
|
| 12 | + |
12 | 13 | import matplotlib.pyplot as plt
|
13 | 14 | import numpy as np
|
14 | 15 | import matplotlib.mlab as mlab
|
15 | 16 | import matplotlib.gridspec as gridspec
|
16 | 17 |
|
17 |
| -# Fixing random state for reproducibility |
18 |
| -np.random.seed(19680801) |
| 18 | +np.random.seed(19680801) # Fix random state for reproducibility. |
19 | 19 |
|
20 | 20 | dt = 0.01
|
21 | 21 | t = np.arange(0, 10, dt)
|
|
30 | 30 | ax0.plot(t, s)
|
31 | 31 | ax1.psd(s, 512, 1 / dt)
|
32 | 32 |
|
33 |
| -plt.show() |
34 |
| - |
35 | 33 | ###############################################################################
|
36 | 34 | # Compare this with the equivalent Matlab code to accomplish the same thing::
|
37 | 35 | #
|
|
57 | 55 | y = 10. * np.sin(2 * np.pi * 4 * t) + 5. * np.sin(2 * np.pi * 4.25 * t)
|
58 | 56 | y = y + np.random.randn(*t.shape)
|
59 | 57 |
|
60 |
| -# Plot the raw time series |
| 58 | +# Plot the raw time series. |
61 | 59 | fig = plt.figure(constrained_layout=True)
|
62 | 60 | gs = gridspec.GridSpec(2, 3, figure=fig)
|
63 | 61 | ax = fig.add_subplot(gs[0, :])
|
64 | 62 | ax.plot(t, y)
|
65 |
| -ax.set_xlabel('time [s]') |
66 |
| -ax.set_ylabel('signal') |
| 63 | +ax.set(xlabel='time [s]', ylabel='signal') |
67 | 64 |
|
68 | 65 | # Plot the PSD with different amounts of zero padding. This uses the entire
|
69 |
| -# time series at once |
| 66 | +# time series at once. |
70 | 67 | ax2 = fig.add_subplot(gs[1, 0])
|
71 | 68 | ax2.psd(y, NFFT=len(t), pad_to=len(t), Fs=fs)
|
72 | 69 | ax2.psd(y, NFFT=len(t), pad_to=len(t) * 2, Fs=fs)
|
73 | 70 | ax2.psd(y, NFFT=len(t), pad_to=len(t) * 4, Fs=fs)
|
74 |
| -ax2.set_title('zero padding') |
| 71 | +ax2.set(title='zero padding') |
75 | 72 |
|
76 |
| -# Plot the PSD with different block sizes, Zero pad to the length of the |
| 73 | +# Plot the PSD with different block sizes, zero pad to the length of the |
77 | 74 | # original data sequence.
|
78 | 75 | ax3 = fig.add_subplot(gs[1, 1], sharex=ax2, sharey=ax2)
|
79 | 76 | ax3.psd(y, NFFT=len(t), pad_to=len(t), Fs=fs)
|
80 | 77 | ax3.psd(y, NFFT=len(t) // 2, pad_to=len(t), Fs=fs)
|
81 | 78 | ax3.psd(y, NFFT=len(t) // 4, pad_to=len(t), Fs=fs)
|
82 |
| -ax3.set_ylabel('') |
83 |
| -ax3.set_title('block size') |
| 79 | +ax3.set(ylabel='', title='block size') |
84 | 80 |
|
85 |
| -# Plot the PSD with different amounts of overlap between blocks |
| 81 | +# Plot the PSD with different amounts of overlap between blocks. |
86 | 82 | ax4 = fig.add_subplot(gs[1, 2], sharex=ax2, sharey=ax2)
|
87 |
| -ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t), noverlap=0, Fs=fs) |
88 |
| -ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t), |
89 |
| - noverlap=int(0.05 * len(t) / 2.), Fs=fs) |
90 |
| -ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t), |
91 |
| - noverlap=int(0.2 * len(t) / 2.), Fs=fs) |
92 |
| -ax4.set_ylabel('') |
93 |
| -ax4.set_title('overlap') |
94 |
| - |
95 |
| -plt.show() |
96 |
| - |
| 83 | +ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t), Fs=fs, noverlap=0) |
| 84 | +ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t), Fs=fs, noverlap=int(0.025*len(t))) |
| 85 | +ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t), Fs=fs, noverlap=int(0.1*len(t))) |
| 86 | +ax4.set(ylabel='', title='overlap') |
97 | 87 |
|
98 | 88 | ###############################################################################
|
99 | 89 | # This is a ported version of a MATLAB example from the signal
|
|
115 | 105 |
|
116 | 106 | ax0.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024,
|
117 | 107 | scale_by_freq=True)
|
118 |
| -ax0.set_title('Periodogram') |
119 |
| -ax0.set_yticks(yticks) |
120 |
| -ax0.set_xticks(xticks) |
| 108 | +ax0.set(title='Periodogram', xticks=xticks, yticks=yticks, ylim=yrange) |
121 | 109 | ax0.grid(True)
|
122 |
| -ax0.set_ylim(yrange) |
123 | 110 |
|
124 | 111 | ax1.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, pad_to=512, noverlap=75,
|
125 | 112 | scale_by_freq=True)
|
126 |
| -ax1.set_title('Welch') |
127 |
| -ax1.set_xticks(xticks) |
128 |
| -ax1.set_yticks(yticks) |
129 |
| -ax1.set_ylabel('') # overwrite the y-label added by `psd` |
| 113 | +ax1.set(title='Welch', xticks=xticks, yticks=yticks, ylim=yrange, |
| 114 | + ylabel='') # overwrite the y-label added by `psd` |
130 | 115 | ax1.grid(True)
|
131 |
| -ax1.set_ylim(yrange) |
132 |
| - |
133 |
| -plt.show() |
134 | 116 |
|
135 | 117 | ###############################################################################
|
136 | 118 | # This is a ported version of a MATLAB example from the signal
|
|
139 | 121 | #
|
140 | 122 | # It uses a complex signal so we can see that complex PSD's work properly.
|
141 | 123 |
|
142 |
| -prng = np.random.RandomState(19680801) # to ensure reproducibility |
143 |
| - |
144 | 124 | fs = 1000
|
145 | 125 | t = np.linspace(0, 0.3, 301)
|
146 | 126 | A = np.array([2, 8]).reshape(-1, 1)
|
147 | 127 | f = np.array([150, 140]).reshape(-1, 1)
|
148 |
| -xn = (A * np.exp(2j * np.pi * f * t)).sum(axis=0) + 5 * prng.randn(*t.shape) |
| 128 | +xn = (A * np.exp(2j * np.pi * f * t)).sum(0) + 5 * np.random.randn(*t.shape) |
149 | 129 |
|
150 | 130 | fig, (ax0, ax1) = plt.subplots(ncols=2, constrained_layout=True)
|
151 | 131 |
|
|
155 | 135 |
|
156 | 136 | ax0.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024,
|
157 | 137 | scale_by_freq=True)
|
158 |
| -ax0.set_title('Periodogram') |
159 |
| -ax0.set_yticks(yticks) |
160 |
| -ax0.set_xticks(xticks) |
| 138 | +ax0.set(title='Periodogram', xticks=xticks, yticks=yticks, ylim=yrange) |
161 | 139 | ax0.grid(True)
|
162 |
| -ax0.set_ylim(yrange) |
163 | 140 |
|
164 | 141 | ax1.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, pad_to=512, noverlap=75,
|
165 | 142 | scale_by_freq=True)
|
166 |
| -ax1.set_title('Welch') |
167 |
| -ax1.set_xticks(xticks) |
168 |
| -ax1.set_yticks(yticks) |
169 |
| -ax1.set_ylabel('') # overwrite the y-label added by `psd` |
| 143 | +ax1.set(title='Welch', xticks=xticks, yticks=yticks, ylim=yrange, |
| 144 | + ylabel='') # overwrite the y-label added by `psd` |
170 | 145 | ax1.grid(True)
|
171 |
| -ax1.set_ylim(yrange) |
172 | 146 |
|
173 | 147 | plt.show()
|
0 commit comments