Skip to content

Commit 2c5ce6e

Browse files
committed
Cleanup psd example.
- Show all figures at once. - Group axes setter calls to have them take less room relative to the psd calls, which are the object of the example. - Align the noverlap examples to make the parallelism between them clearer. - Don't introduce a second random state with a single use.
1 parent 165df1b commit 2c5ce6e

File tree

1 file changed

+21
-47
lines changed

1 file changed

+21
-47
lines changed
+21-47
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
"""
22
========
3-
Psd Demo
3+
PSD Demo
44
========
55
66
Plotting Power Spectral Density (PSD) in Matplotlib.
@@ -9,13 +9,13 @@
99
many useful libraries for computing a PSD. Below we demo a few examples
1010
of how this can be accomplished and visualized with Matplotlib.
1111
"""
12+
1213
import matplotlib.pyplot as plt
1314
import numpy as np
1415
import matplotlib.mlab as mlab
1516
import matplotlib.gridspec as gridspec
1617

17-
# Fixing random state for reproducibility
18-
np.random.seed(19680801)
18+
np.random.seed(19680801) # Fix random state for reproducibility.
1919

2020
dt = 0.01
2121
t = np.arange(0, 10, dt)
@@ -30,8 +30,6 @@
3030
ax0.plot(t, s)
3131
ax1.psd(s, 512, 1 / dt)
3232

33-
plt.show()
34-
3533
###############################################################################
3634
# Compare this with the equivalent Matlab code to accomplish the same thing::
3735
#
@@ -57,43 +55,35 @@
5755
y = 10. * np.sin(2 * np.pi * 4 * t) + 5. * np.sin(2 * np.pi * 4.25 * t)
5856
y = y + np.random.randn(*t.shape)
5957

60-
# Plot the raw time series
58+
# Plot the raw time series.
6159
fig = plt.figure(constrained_layout=True)
6260
gs = gridspec.GridSpec(2, 3, figure=fig)
6361
ax = fig.add_subplot(gs[0, :])
6462
ax.plot(t, y)
65-
ax.set_xlabel('time [s]')
66-
ax.set_ylabel('signal')
63+
ax.set(xlabel='time [s]', ylabel='signal')
6764

6865
# Plot the PSD with different amounts of zero padding. This uses the entire
69-
# time series at once
66+
# time series at once.
7067
ax2 = fig.add_subplot(gs[1, 0])
7168
ax2.psd(y, NFFT=len(t), pad_to=len(t), Fs=fs)
7269
ax2.psd(y, NFFT=len(t), pad_to=len(t) * 2, Fs=fs)
7370
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')
7572

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
7774
# original data sequence.
7875
ax3 = fig.add_subplot(gs[1, 1], sharex=ax2, sharey=ax2)
7976
ax3.psd(y, NFFT=len(t), pad_to=len(t), Fs=fs)
8077
ax3.psd(y, NFFT=len(t) // 2, pad_to=len(t), Fs=fs)
8178
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')
8480

85-
# Plot the PSD with different amounts of overlap between blocks
81+
# Plot the PSD with different amounts of overlap between blocks.
8682
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')
9787

9888
###############################################################################
9989
# This is a ported version of a MATLAB example from the signal
@@ -115,22 +105,14 @@
115105

116106
ax0.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024,
117107
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)
121109
ax0.grid(True)
122-
ax0.set_ylim(yrange)
123110

124111
ax1.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, pad_to=512, noverlap=75,
125112
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`
130115
ax1.grid(True)
131-
ax1.set_ylim(yrange)
132-
133-
plt.show()
134116

135117
###############################################################################
136118
# This is a ported version of a MATLAB example from the signal
@@ -139,13 +121,11 @@
139121
#
140122
# It uses a complex signal so we can see that complex PSD's work properly.
141123

142-
prng = np.random.RandomState(19680801) # to ensure reproducibility
143-
144124
fs = 1000
145125
t = np.linspace(0, 0.3, 301)
146126
A = np.array([2, 8]).reshape(-1, 1)
147127
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)
149129

150130
fig, (ax0, ax1) = plt.subplots(ncols=2, constrained_layout=True)
151131

@@ -155,19 +135,13 @@
155135

156136
ax0.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024,
157137
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)
161139
ax0.grid(True)
162-
ax0.set_ylim(yrange)
163140

164141
ax1.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, pad_to=512, noverlap=75,
165142
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`
170145
ax1.grid(True)
171-
ax1.set_ylim(yrange)
172146

173147
plt.show()

0 commit comments

Comments
 (0)