|
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 | 16 | importmatplotlib.gridspecasgridspec |
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() |