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