直接上代码，理论可以去知乎看。
#Import necessary libraries
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn as snsimport pywt
from scipy.ndimage import gaussian_filter1d
from scipy.signal import chirp
import matplotlib.gridspec as gridspec
from scipy import signal
from skimage import filters,img_as_float
from skimage.io import imread, imshow
from skimage.color import rgb2hsv, rgb2gray, rgb2yuv
from skimage import color, exposure, transform
from skimage.exposure import equalize_hist
from scipy import fftpack, ndimageHow to choose scale for wavelet transform
t_min=0
t_max=10
fs=100
dt = 1/fs
time = np.linspace(t_min, t_max, 1500)
#To understand the behaviour of scale, we used a smooth constant signal with a discontinuity. Adding discontinuity to the constant will have a rectangular shape.
w = chirp(time, f0=10, f1=50, t1=10, method='quadratic')#Compute Wavelet Transform
scale = [10,20,30,50,100]#Plot signal, FFT, and scalogram(to represent wavelet transform)
fig,axes = plt.subplots(nrows=1,ncols=5,figsize=(25,4))
for i in range(2):for j in range(5):#Scalogramscales = np.arange(1,scale[j],1)coef,freqs = pywt.cwt(w,scales,'morl')freqs = pywt.scale2frequency('morl',scales,precision=8)if i == 0:axes[j].set_title("Scalogram from scale {} to {}".format(1,scale[j]))if i == 0:axes[j].pcolormesh(time, scales, coef,cmap='Greys')axes[j].set_ylabel("Scale")
plt.show();
scales = np.arange(1,20,1)
coef,freqs = pywt.cwt(w,scales,'morl',1/fs)
fig,axes = plt.subplots(nrows=1,ncols=2,figsize=(12,5))
axes[0].set_title("Scalogram")
axes[0].pcolormesh(time, scales, coef,cmap='Greys')
axes[0].set_xlabel("Time")
axes[0].set_ylabel("Scale")
axes[1].set_title("Spectrogram")
axes[1].pcolormesh(time, freqs, coef,cmap='Greys')
axes[1].set_xlabel("Time")
axes[1].set_ylabel("Pseudo Frequency")
plt.show();
families = ['gaus1','gaus2','gaus3','gaus4','gaus5','gaus6','gaus7','gaus8','mexh','morl']
cols = 5
rows = 4
scales = np.arange(1,20,1)
fig,axes = plt.subplots(nrows = rows,ncols=5,figsize=(3*cols,2*rows))
fig.tight_layout(pad=1.0, w_pad=1.0, h_pad=3)
for i,family in enumerate(families):c = i%5r = round(i//5)coef,freqs = pywt.cwt(w,scales,family,1/fs)psi, x = pywt.ContinuousWavelet(family).wavefun(level=10)axes[r*2,c].set_title(family)axes[(r*2)+1,c].pcolormesh(time, freqs, coef,cmap='Blues')axes[(r*2)+1,c].set_xlabel("Time")axes[(r*2)+1,c].set_ylabel("Scale")axes[r*2,c].plot(x, psi)axes[r*2,c].set_xlabel("X")axes[r*2,c].set_ylabel("Psi")
Constant Signal
fs = 100 #Sampling frequency
time = np.arange(-3,3,1/fs) #create time
n = len(time)
T=1/fs
print("We consider {} samples".format(n))
constant = np.ones(n) #Amblitude will be one(constant value)
freq = np.linspace(-1.0/(2.0*T), 1.0/(2.0*T), n)#Compute Fourier transform of Constant signal
fft = fftpack.fft(constant)
freq = fftpack.fftfreq(time.shape[0],T)
phase = np.angle(fft)
phase = phase / np.pi#Compute Wavelet Transform
scales = np.arange(1,6,1)
coef,freqs = pywt.cwt(constant,scales,'gaus1')#Plot signal, FFT, and scalogram(to represent wavelet transform)
fig,axes = plt.subplots(ncols=3,figsize=(18,4))#Signal
axes[0].set_title("Constant")
axes[0].plot(time, constant)
axes[0].set_xlabel("Time")
axes[0].set_ylabel("Amplitude")#Fourier
axes[1].set_title("Fourier Transform")
axes[1].plot(freq, np.abs(fft)/n)
axes[1].set_xlabel("Frequency")
axes[1].set_ylabel("Magnitude")#Scalogram
axes[2].set_title("Scalogram")
axes[2].pcolormesh(time, scales, coef,cmap='bone')
axes[2].set_xlabel("Time")
axes[2].set_ylabel("Scale")
plt.show();
Adding discontinuity to constant signal
constant[300:340]=0#Compute Fourier transform of Constant signal
fft = fftpack.fft(constant)
phase = np.angle(fft)
phase = phase / np.pi#Compute Wavelet Transform
scales = np.arange(1,6,1)
coef,freqs = pywt.cwt(constant,scales,'gaus1')#Plot signal, FFT, and scalogram(to represent wavelet transform)
fig,axes = plt.subplots(ncols=3,figsize=(18,4))#Signal
axes[0].set_title("Constant")
axes[0].plot(time, constant)
axes[0].set_xlabel("Time")
axes[0].set_ylabel("Amplitude")#Fourier
axes[1].set_title("Fourier Transform")
axes[1].plot(freq, np.abs(fft)/n)
axes[1].set_xlabel("Frequency")
axes[1].set_ylabel("Magnitude")#Scalogram
axes[2].set_title("Scalogram")
axes[2].pcolormesh(time, scales, coef,cmap='bone')
axes[2].set_xlabel("Time")
axes[2].set_ylabel("Scale")
plt.show();
Rectangular Pulse
N = 50000 #number of samples
fs = 1000 #sample frequency
T = 1/fs #interval
time = np.linspace(-(N*T), N*T, N)
rect = np.zeros(time.shape)
for i in range(time.shape[0]):if time[i] > -0.5 and time[i] < 0.5:rect[i] = 1.0
print("We consider {} samples".format(N))
freq = np.linspace(-1.0/(2.0*T), 1.0/(2.0*T), N)#compute Fourier Trainsform
fft_rect = np.fft.fft(rect)
fr = np.fft.fftfreq(N)
phase = np.angle(fft_rect)
phase = phase / np.pi
freqrect = np.fft.fftfreq(time.shape[-1])
fft_rect = np.fft.fftshift(fft_rect)#compute wavelet transform
scales = np.arange(1,25,1)
coef,freqs = pywt.cwt(rect,scales,'gaus1')#Plot
#signal
fig,axes = plt.subplots(ncols=3,figsize=(21,5))
axes[0].set_title("Rectangular signal")
axes[0].plot(time, rect)
axes[0].set_xlim(-1,1)
axes[0].set_xlabel("Time")
axes[0].set_ylabel("rectangular pulse")#Fourier transform
axes[1].set_title("Fourier Transform")
axes[1].plot(freq,np.abs(fft_rect)*2/fs)
axes[1].set_xlim(-40,40)
axes[1].set_xlabel("Frequency")
axes[1].set_ylabel("Magnitude")#wavelet
axes[2].set_title("Scalogram ")
axes[2].pcolormesh(time, scales, coef,cmap='bone')
axes[2].set_xlim(-2,2)
axes[2].set_xlabel("Time")
axes[2].set_ylabel("Scale")
plt.show();
Sine and cosine waves
fs = 1000 #sampling frequency
interval = 1/fs #sampling interval
t_min = -1 #start time
t_max = 1 # end time
dt=1/fs
time = np.arange(t_min,t_max,interval)n = len(time)
print("We consider {} samples".format(n))f = (fs/2)*np.linspace(0,1,int(n/2)) #frequencyfreq = [200,130] #signal frequencies
scales1 = np.arange(1,20,1)#Create signal with 200 hz frequency
sinewave1 = np.sin(2*np.pi*freq[0]*time)
new = sinewave1/np.square(time)
#compute fourier transform
fft1 = np.fft.fft(sinewave1)
fr = np.fft.fftfreq(n, d=dt)
phase = np.angle(fft1)
phase = phase / np.pi
fft1 = fft1[0:int(n/2)]#compute wavelet
coef1,freqs1 = pywt.cwt(sinewave1,scales1,'morl')#plot
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)
ax = plt.subplot(gs[0, :])
ax.set_title("Sinusoidal Signal - 200 Hz")
ax.plot(time,sinewave1)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax2 = plt.subplot(gs[1, 0])
ax2.plot(f,np.abs(fft1)*2/fs)
ax2.set_xlabel("Frequency")
ax2.set_ylabel("DFT values")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs1/dt, coef1)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")
plt.show;
scales = np.arange(1,20,1)
#Create signal with 130 hz frequency
sinewave2 = np.sin(2*np.pi*freq[1]*time)#compute fourier transform
fft2 = np.fft.fft(sinewave2)
fft2=fft2[0:int(n/2)]#compute wavelet
coef2,freqs2 = pywt.cwt(sinewave2,scales,'morl')
#plot
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)
ax = plt.subplot(gs[0, :])
ax.set_title("Sinusoidal Signal - 130 Hz")
ax.plot(time,sinewave2)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax2 = plt.subplot(gs[1, 0])
ax2.plot(f,np.abs(fft2)*2/fs)
ax2.set_xlim(0,300)
ax2.set_xlabel("Frequency")
ax2.set_ylabel("Magnitude")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs2/dt, coef2)
ax3.set_xlabel("Time")
ax3.set_ylabel("Scale")
plt.show();
scales = np.arange(1,30,1)sum = sinewave1+sinewave2
fft3 = np.fft.fft(sum)
fft3=fft3[0:int(n/2)]
coef3,freqs3 = pywt.cwt(sum,scales,'morl')#plot
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)ax = plt.subplot(gs[0, :])
ax.set_title("Sum of Sinusoidal")
ax.plot(time,sum)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax2 = plt.subplot(gs[1, 0])
ax2.plot(f,np.abs(fft3)*2/fs)
ax2.set_xlabel("Frequency")
ax2.set_ylabel("Magnitude")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs3/dt, coef3)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")
plt.show();
Non stationary signals
size = len(time)//3
scales = np.arange(1,31,1)
sig = np.zeros(time.shape)
sig[:size]=np.sin(2*np.pi*200*time[:size])
sig[size:size*2]=np.sin(2*np.pi*130*time[size:size*2])
sig[size*2:]=np.cos(2*np.pi*50*time[size*2:])
fft = np.fft.fft(sig)
fft=fft[0:int(n/2)]
coef,freqs = pywt.cwt(sig,scales,'gaus8')
stft_f, stft_t, Sxx = signal.spectrogram(sig, fs,window='hann', nperseg=64)
#plot
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)ax = plt.subplot(gs[0, 0])
ax.set_title("Sinusoidal Signal- Frequency vary over time")
ax.plot(time,sig)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax1 = plt.subplot(gs[0, 1])
ax1.plot(f,np.abs(fft)*2/fs)
ax1.set_xlabel("Frequency")
ax1.set_ylabel("Magnitude")
ax2 = plt.subplot(gs[1, 0])
ax2.pcolormesh(stft_t, stft_f, Sxx)
ax2.set_ylabel("Frequency")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs/dt, coef)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")
Linear Chirp Signal
def plot_chirp_transforms(type_,f0,f1):#Create linear chirp signal with frequency between 50Hz and 10Hzt_min=0t_max=10time = np.linspace(t_min, t_max, 1500)N = len(time)interval = (t_min+t_max)/Nfs = int(1/interval)dt=1/fsf = (fs/2)*np.linspace(0,1,int(N/2))w1 = chirp(time, f0=f0, f1=f1, t1=10, method=type_.lower())#Compute FFTw1fft = np.fft.fft(w1)w1fft=w1fft[0:int(N/2)]#Compute Wavelet transformscales=np.arange(1,50,1)wcoef,wfreqs = pywt.cwt(w1,scales,'morl')#Compute Short Time Fourier transfomrstft_f, stft_t, Sxx = signal.spectrogram(w1, fs,window='hann', nperseg=64,noverlap=32)#Plot the resultsgs = gridspec.GridSpec(2,2)gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)ax = plt.subplot(gs[0, 0])ax.set_title("Chirp - "+type_+" ({}Hz to {}Hz)".format(f0,f1))ax.plot(time,w1)ax.set_xlabel("Time(s)")ax.set_ylabel("Amplitude")ax1 = plt.subplot(gs[0, 1])ax1.plot(f,np.abs(w1fft)*2/fs)plt.grid()ax1.set_title("FFT - "+type_+" chirp signal ({}Hz to {}Hz)".format(f0,f1))ax1.set_xlabel("Frequency")ax1.set_ylabel("Magnitude")ax2 = plt.subplot(gs[1, 0])ax2.set_title("STFT - "+type_+" chirp signal ({}Hz to {}Hz)".format(f0,f1))ax2.pcolor(stft_t, stft_f, Sxx,cmap='copper')ax2.set_xlabel("Time")ax2.set_ylabel("Frequency")ax3 = plt.subplot(gs[1, 1])ax3.set_title("WT - "+type_+" chirp signal ({}Hz to {}Hz)".format(f0,f1))ax3.pcolor(time, wfreqs/dt, wcoef,cmap='copper')ax3.set_ylim(5,75)ax3.set_xlabel("Time")ax3.set_ylabel("Frequency")
plot_chirp_transforms('Linear',50,10)
plot_chirp_transforms('Linear',10,50)
plot_chirp_transforms('Quadratic',50,10)
Trapezoid
N = 5000 #number of samples
fs = 1000 #sample frequency
T = 1/fs #interval
time = np.linspace(-5, 5, N)
trapzoid_signal = (time*np.where(time>0,1,0))-((time-1)*np.where((time-1)>0,1,0))-((time-2)*np.where((time-2)>0,1,0))+((time-3)*np.where((time-3)>0,1,0))#tra = trapzoid_signal(time)scales = np.arange(1,51,1)
coef,freqs = pywt.cwt(trapzoid_signal,scales,'gaus1')
#compute Fourier Trainsform
fft = np.fft.fft(trapzoid_signal)
freq = np.fft.fftfreq(time.shape[-1],T)
fftShift = np.fft.fftshift(fft)
freqShift=np.fft.fftshift(freq)#Plot signal and FFT
fig,axes = plt.subplots(nrows=2,ncols=3,figsize=(24,10))
axes[0,0].set_title("Trapezoidal signal")
axes[0,0].plot(time, trapzoid_signal)
axes[0,0].set_xlabel("Time")
axes[0,0].set_ylabel("Trapezoidal pulse")
axes[0,1].set_title("Fourier Transform - trapezoidal")
axes[0,1].plot(freqShift,np.abs(fftShift)*2/fs)
axes[0,1].set_xlim(-20,20)
axes[0,1].set_xlabel("Frequency")
axes[0,1].set_ylabel("Magnitude")
axes[0,2].set_title("Scalogram - trapezoidal")
axes[0,2].pcolor(time,scales,coef,cmap='BrBG')
axes[0,2].set_xlim(-5,5)
axes[0,2].set_xlabel("Time")
axes[0,2].set_ylabel("Scale")trapzoid_signal = ((time+1)*np.where((time+1)>0,1,0))-(time*np.where(time>0,1,0))-((time-1)*np.where((time-1)>0,1,0))+((time-2)*np.where((time-2)>0,1,0))
coef,freqs = pywt.cwt(trapzoid_signal,scales,'gaus1')
#compute Fourier Trainsform
fft = np.fft.fft(trapzoid_signal)
freq = np.fft.fftfreq(time.shape[-1],T)
fftShift = np.fft.fftshift(fft)
freqShift=np.fft.fftshift(freq)#Plot signal and FFT
axes[1,0].plot(time, trapzoid_signal)
axes[1,0].set_xlabel("Time")
axes[1,0].set_ylabel("Trapezoidal pulse")
axes[1,1].plot(freqShift,np.abs(fftShift)*2/fs)
axes[1,1].set_xlim(-20,20)
axes[1,1].set_xlabel("Frequency")
axes[1,1].set_ylabel("Magnitude")
axes[1,2].pcolor(time,scales,coef,cmap='BrBG')
axes[1,2].set_xlim(-5,5)
axes[1,2].set_xlabel("Time")
axes[1,2].set_ylabel("Scale")
#Image Generation
sigma = 20
fake_image = np.repeat(a=np.repeat(a=160,repeats=512),repeats=512).reshape([512,512])
fake_image_translated = fake_image.copy()
fake_image[250:350,200:300] = 190
fake_image_translated[250:350,250:350] = 190
fake_image = filters.gaussian(fake_image, sigma=sigma, preserve_range=True)
fake_image_translated = filters.gaussian(fake_image_translated, sigma=sigma, preserve_range=True)
fake_image = rgb2gray(fake_image)/255
fake_image_translated = rgb2gray(fake_image_translated)/255#Fourier transform
fft_fake = fftpack.fft2(fake_image)
fft_fake = fftpack.fftshift(fft_fake)
fft_fake2 = fftpack.fft2(fake_image_translated)
fft_fake2 = fftpack.fftshift(fft_fake2)#Plot
fig,axes = plt.subplots(ncols=4,figsize=(16,4));
axes[0].set_title("Fake Image");
axes[0].imshow(fake_image,cmap='gray');
axes[1].set_title("FFT");
axes[1].imshow(np.log(np.abs(fft_fake)),cmap='gray');
axes[2].set_title("Shifted Image");
axes[2].imshow(fake_image_translated,cmap='gray');
axes[3].set_title("FFT-shifted image");
axes[3].imshow(np.log(np.abs(fft_fake2)),cmap='gray');
plt.show();
titles = ['Approximation', ' Horizontal detail','Vertical detail', 'Diagonal detail']
coeffs_fi = pywt.dwt2(fake_image, 'haar')
coeffs_fi_translated = pywt.dwt2(fake_image_translated, 'haar')
coef_array = []
cA1, (cH1, cV1, cD1) = coeffs_fi
cA2, (cH2, cV2, cD2) = coeffs_fi_translated
coef_array.append([cA1,cH1, cV1, cD1])
coef_array.append([cA2,cH2, cV2, cD2])
fig,axes = plt.subplots(nrows=2,ncols=5,figsize=(25,10))
axes[0,0].set_title("Image")
axes[0,0].imshow(fake_image,cmap='gray')
axes[1,0].set_title("Image")
axes[1,0].imshow(fake_image_translated,cmap='gray')
for i,arr in enumerate(coef_array):for idx,coef in enumerate(arr):axes[i,idx+1].set_title(titles[idx])axes[i,idx+1].imshow(coef,cmap='gray')
from scipy import ndimage#Image Generation
line_image = np.repeat(a=np.repeat(a=0,repeats=512),repeats=512).reshape([512,512])
line_image_translated = line_image.copy()
line_image[192:320,128:354] = 165
line_image[224:288,192:320] = 70
line_image_translated[128:354,192:320] = 165
line_image_translated[192:320,224:288] = 70
line_image = filters.gaussian(line_image, sigma=sigma, preserve_range=True)
line_image_translated = filters.gaussian(line_image_translated, sigma=sigma, preserve_range=True)
line_image = rgb2gray(line_image)/255
line_image_translated = rgb2gray(line_image_translated)/255
line_image_rot = ndimage.rotate(line_image, 45, reshape=False)#Fourier transform
fft_line = fftpack.fft2(line_image)
fft_line = fftpack.fftshift(fft_line)
fft_line_t = fftpack.fft2(line_image_translated)
fft_line_t = fftpack.fftshift(fft_line_t)
fft_line_r = fftpack.fft2(line_image_rot)
fft_line_r = fftpack.fftshift(fft_line_r)#Plot
fig,axes = plt.subplots(ncols=6,figsize=(24,4))
axes[0].set_title("Fake Image")
axes[0].imshow(line_image,cmap='gray')
axes[1].set_title("FFT")
axes[1].imshow(np.log(np.abs(fft_line)),cmap='gray')
axes[2].set_title("Shifted Image")
axes[2].imshow(line_image_translated,cmap='gray')
axes[3].set_title("FFT-shifted image")
axes[3].imshow(np.log(np.abs(fft_line_t)),cmap='gray')
axes[4].set_title("Rotated image")
axes[4].imshow(line_image_rot,cmap='gray')
axes[5].set_title("FFT-Rotated")
axes[5].imshow(np.log(np.abs(fft_line_r)),cmap='gray')
plt.show();
coef_line = pywt.dwt2(line_image, 'haar')
coef_line_t = pywt.dwt2(line_image_translated, 'haar')
coef_line_r = pywt.dwt2(line_image_rot, 'haar')
dwt_coef_array = []
lcA1, (lcH1, lcV1, lcD1) = coef_line
lcA2, (lcH2, lcV2, lcD2) = coef_line_t
lcA3, (lcH3, lcV3, lcD3) = coef_line_r
dwt_coef_array.append([lcA1,lcH1, lcV1, lcD1])
dwt_coef_array.append([lcA2,lcH2, lcV2, lcD2])
dwt_coef_array.append([lcA3,lcH3, lcV3, lcD3])
fig,axes = plt.subplots(nrows=3,ncols=5,figsize=(25,15))
axes[0,0].set_title("Fake Image")
axes[0,0].imshow(line_image,cmap='gray')
axes[1,0].set_title("Translated Image")
axes[1,0].imshow(line_image_translated,cmap='gray')
axes[2,0].set_title("Rotated Image")
axes[2,0].imshow(line_image_rot,cmap='gray')
for i,arr in enumerate(dwt_coef_array):for idx,coef in enumerate(arr):axes[i,idx+1].set_title(titles[idx])axes[i,idx+1].imshow(coef,cmap='gray')
from PIL import Image
# open the original image
original_img = Image.open("/content/drive/MyDrive/DSIP/parrot1.jpg")#rotate image
rot_180 = original_img.rotate(180, Image.NEAREST, expand = 1)# close all our files objectI = np.array(original_img)
I_rot = np.array(rot_180)original_img.close()I_grey = rgb2gray(I)
I_rot_grey = rgb2gray(I_rot)fft2 = fftpack.fft2(I_grey)
fftshift = fftpack.fftshift(fft2)
fftrot2 = fftpack.fft2(I_rot_grey)
fftrotshift = fftpack.fftshift(fftrot2)coeffs2 = pywt.dwt2(I_grey, 'haar')
cA, (cH, cV, cD) = coeffs2
titles = ['Approximation', ' Horizontal detail','Vertical detail', 'Diagonal detail']
coeffs3 = pywt.dwt2(I_rot_grey, 'haar')
cA1, (cH1, cV1, cD1) = coeffs3
fig,axes = plt.subplots(ncols=6,nrows=2,figsize=(24,8))axes[0,0].set_title("Image")
axes[0,0].imshow(img_as_float(I_grey),cmap='gray')
axes[0,1].set_title("FFT")
axes[0,1].imshow(np.log(np.abs(fftshift)),cmap='gray')
axes[1,0].set_title("Flip Image")
axes[1,0].imshow(img_as_float(I_rot_grey),cmap='gray')
axes[1,1].set_title("FFT-flip image")
axes[1,1].imshow(np.log(np.abs(fftrotshift)),cmap='gray')for idx,coef in enumerate((cA,cH,cV,cD)):axes[0,idx+2].set_title(titles[idx])axes[0,idx+2].imshow(coef,cmap='gray')
for idx,coef in enumerate((cA1,cH1,cV1,cD1)):axes[1,idx+2].set_title(titles[idx])axes[1,idx+2].imshow(coef,cmap='gray')
plt.show();
知乎学术咨询：
https://www.zhihu.com/consult/people/792359672131756032?isMe=1
擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。
