有趣的傅里叶变换与小波变换对比(Python)

不严谨的说,时域和频域分析就是在不同的空间看待问题的,不同空间所对应的原子(基函数)是不同的。你想一下时域空间的基函数是什么?频域空间的基函数是什么?一般的时-频联合域空间的基函数是什么?小波域空间的基函数是什么?

有的空间域比较容易分析,有的空间域不容易分析。

举个例子吧,首先加载一个双曲Chirp信号,数据的采样频率为2048Hz,第一个Chirp信号持续时间为0.1~0.68秒,第二个Chirp信号持续时间为0.1~0.75 秒,第一个Chirp信号在时间t处的瞬时频率为(单位Hz):

第二个Chirp信号在时间t处的瞬时频率为(单位Hz):

看一下从时域空间看待的时域图

然后看一下频域空间的频谱图
傅里叶变换(FT)比较擅长识别信号中存在的频率分量, 但是FT无法定位频率分量。绘制上面信号的幅值谱,并放大0到200Hz之间的区域

再看一下一般的时频域空间的时频谱图,以短时傅里叶变换为例
傅里叶变换不提供时间信息,为了定位频率,短时傅里叶变换STFT方法将信号分割成不同的窗,并对每个窗执行FT。STFT的时频分析窗口如下:

STFT提供了信号时间-频率域中的一些信息, 但是选择窗的大小是关键。对于STFT时频分析,选择更短的窗以牺牲频率分辨率为代价从而获得良好的时间分辨率。相反,选择较大的窗以时间分辨率为代价从而获得良好的频率分辨率(著名的测不准原理)。一旦STFT的分析窗确定后,将在整个分析中保持不变(最致命的缺陷)。以 200 毫秒的时间窗大小绘制上述双曲Chirp信号的频谱图,频谱图上的瞬时频率为黑色虚线段。

然后绘制时间窗大小为50毫秒的频谱图

两个图的结果是显而易见的,没有单一的窗口大小可以解析此类信号的整个频率信息。
最后看一下小波空间对应的小波时频谱图
连续小波变换 CWT是为了克服 STFT中固有的时频分辨率问题。CWT的时频分析窗口如下:

CWT和人类的听觉系统非常一致:在低频处有更好的频率定位能力,在高频处有更好的时间定位能力。绘制 CWT时尺度谱(尺度谱是作为时间和频率绘制的 CWT的绝对值),因为CWT 中的频率是对数的,所以使用对数频率轴。

从图中可以清楚地看出信号中两个双曲Chirp信号的存在,CWT可以比较准确估计持续时间的瞬时频率,而无需担心选择窗的大小。要了解小波系数幅度增长速度有多快,可以看一下3-D 图

在尺度谱上绘制一下瞬时频率,可见瞬时频率与尺度谱特征非常吻合

看到了吧,从不同的空间域(角度)看待问题,分析的难度也不一样。

开始正题Wavelet vs Fourier transform

#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, ndimage
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")

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();

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();

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();

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();

scales = np.arange(1,40,1)
#creating two cosine waves at 50 hz frequency
coswave = np.cos(2*np.pi*50*time)
#compute Fourier Trainsform of amblitude with 10 Hz
fft = np.fft.fft(coswave)
fft=fft[0:int(n/2)]
coef,freqs = pywt.cwt(coswave,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("Cosine Signal - 50 Hz")
ax.plot(time,coswave)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax2 = plt.subplot(gs[1, 0])
ax2.plot(f,np.abs(fft)*2/fs)
ax2.set_xlabel("Frequency")
ax2.set_ylabel("Magnitude")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs/dt, coef)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")
plt.show()

Inject a Signal to the Sum signal

scales = np.arange(1,40,1)
sum_disc = sum+np.real(np.exp(-50*(time-0.4)**2)*np.exp(1j*2*np.pi*400*(time-0.4)))
fft = np.fft.fft(sum_disc)
fft=fft[0:int(n/2)]
coef,freqs = pywt.cwt(sum_disc,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("Signal with variation")
ax.plot(time,sum_disc)
ax.set_xlabel("Time(s)")
ax.set_ylabel("Amplitude")
ax2 = plt.subplot(gs[1, 0])
ax2.plot(f,np.abs(fft)*2/fs)
ax2.set_xlim(0,500)
ax2.set_xlabel("Frequency")
ax2.set_ylabel("Magnitude")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, freqs*1000, coef)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")

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('Logarithmic',50,10)

plot_chirp_transforms('Hyperbolic',50,10)

Triangular Signal

N = 1000 #number of samples
fs = 1000 #sample frequency
T = 1/fs #interval
time = np.linspace(-2, 2, N)
tri = np.where(np.abs(time)<=.5,1,0)
tri = np.where(tri==1,.5-np.abs(time),0)print("We consider {} samples".format(N))
scales = np.arange(1,51,1)
coef,freqs = pywt.cwt(tri,scales,'gaus1')
#compute Fourier Trainsform
fft = np.fft.fft(tri)
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("Triangular signal")
axes[0,0].plot(time, tri)
axes[0,0].set_xlabel("Time")
axes[0,0].set_ylabel("Triangular pulse")
axes[0,1].set_title("Fourier Transform")
axes[0,1].plot(freqShift,np.abs(fftShift)*2/fs)
axes[0,1].set_xlabel("Frequency")
axes[0,1].set_xlim(-50,50)
axes[0,1].set_ylabel("Magnitude")
axes[0,2].set_title("Scalogram")
axes[0,2].pcolor(time,scales,coef,cmap='copper')
axes[0,2].set_xlim(-2,2)
axes[0,2].set_xlabel("Time")
axes[0,2].set_ylabel("Scale")indices = np.where(tri>0)[0]
new_indices = indices+120
temp = tri.copy()
tri = np.zeros(tri.shape)
tri[new_indices] = temp[indices]coef,freqs = pywt.cwt(tri,scales,'gaus1')
#compute Fourier Trainsform
fft = np.fft.fft(tri)
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, tri)
axes[1,0].set_xlabel("Time")
axes[1,0].set_ylabel("Triangular pulse")
axes[1,1].plot(freqShift,np.abs(fftShift)*2/fs)
axes[1,1].set_xlim(-50,50)
axes[1,1].set_xlabel("Frequency")
axes[1,1].set_ylabel("Magnitude")
axes[1,2].pcolor(time,scales,coef,cmap='copper')
axes[1,2].set_xlim(-2,2)
axes[1,2].set_xlabel("Time")
axes[1,2].set_ylabel("Scale")

Audio

from scipy.io import wavfile
import scipy
#Read Audio and compute time
sr, data = wavfile.read('5.wav')
dt = 1/sr
time = np.arange(0,1,dt)
#Find FFT and frequencies
fft_aud = np.fft.fft(data)
fft_aud=fft_aud[0:int(sr/2)]
freq = (sr/2)*np.linspace(0,1,int(sr/2))
plt.plot(time,data)#Compute STFT
stft_f, stft_t, Sxx = signal.spectrogram(data, sr,window='hann', nperseg=256,noverlap=64)#Compute Wavelet Transform(morlet)
widths = np.arange(1, 31)
wt,wfreqs = pywt.cwt(data,widths,'morl')
gs = gridspec.GridSpec(2,2)
gs.update(left=0, right=4,top=2,bottom=0, hspace=.2,wspace=.1)ax0 = plt.subplot(gs[0, 0])
ax0.set_title("Audio Signal")
ax0.plot(time,data)
ax0.set_xlabel("Time")
ax0.set_ylabel("Amplitude")
ax1 = plt.subplot(gs[0, 1])
ax1.plot(freq, np.abs(fft_aud))
ax1.set_xlabel("Frequency")
ax1.set_ylabel("Magnitude")
ax2 = plt.subplot(gs[1, 0])
ax2.pcolor(stft_t, stft_f, Sxx, cmap='copper')
ax2.set_xlabel("Time")
ax2.set_ylabel("Frequency")
ax3 = plt.subplot(gs[1, 1])
ax3.pcolormesh(time, wfreqs/dt, wt)
ax3.set_xlabel("Time")
ax3.set_ylabel("Frequency")
plt.show()

from PIL import Image
# open the original image
original_img = Image.open("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();

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")
知乎学术咨询:
https://www.zhihu.com/consult/people/792359672131756032?isMe=1
axes[1,2].set_ylabel("Scale")

工学博士,担任《Mechanical System and Signal Processing》《中国电机工程学报》《控制与决策》等期刊审稿专家,擅长领域:现代信号处理,机器学习,深度学习,数字孪生,时间序列分析,设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/351199.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

Arduino入门2——常用函数及用法

Arduino入门2——串口驱动函数及用法 IO串口 上期&#xff0c;我们简单的认识了一下Arduino&#xff0c;浅浅的入了个门&#xff0c;这一期我们介绍以下Arduino串口常用的函数及用法 IO 常用串口库函数如下&#xff1a; 函数名用法及解析pinMode()用于IO口初始化digitalWrite…

28.启动与暂停程序

上一个内容&#xff1a;27.设计注入功能界面 以它 27.设计注入功能界面 的代码为基础进行修改 点击添加游戏按钮之后就把游戏启动了 CWndINJ.cpp文件中修改&#xff1a; void CWndINJ::OnBnClickedButton1() {// TODO: 在此添加控件通知处理程序代码/*ExeLst.InsertItem(0, L…

Opus从入门到精通(五)OggOpus封装器全解析

Opus从入门到精通(五)OggOpus封装器全解析 为什么要封装 前面Opus从入门到精通(四)Opus解码程序实现提到如果不封装会有两个问题: 无法从文件本身获取音频的元数据(采样率,声道数,码率等)缺少帧分隔标识,无法从连续的文件流中分隔帧(尤其是vbr情况) 针对上面的问题我们可以…

4-异常-log4j配置日志滚动覆盖出现日志丢失问题

4-异常-log4j配置日志打印滚动覆盖出现日志丢失问题(附源码分析) 更多内容欢迎关注我&#xff08;持续更新中&#xff0c;欢迎Star✨&#xff09; Github&#xff1a;CodeZeng1998/Java-Developer-Work-Note 技术公众号&#xff1a;CodeZeng1998&#xff08;纯纯技术文&…

springboot依赖管理和自动配置

依赖管理和自动配置 依赖管理和自动配置依赖管理什么是依赖管理修改自动仲裁/默认版本号 starter场景启动器starter场景启动器基本介绍官方提供的starter第三方starter 自动配置自动配置基本介绍SpringBoot自动配置了哪些?如何修改默认配置如何修改默认扫描包结构resources\ap…

openGauss 6.0.0 一主二备集群安装及使用zcbus实现Oracle到openGauss的数据同步

一、前言 openGauss 6.0.0-RC1是openGauss 2024年3月发布的创新版本&#xff0c;该版本生命周期为0.5年。根据openGauss官网介绍&#xff0c;6.0.0-RC1与之前的版本特性功能保持兼容,另外&#xff0c;在和之前版本兼容的基础上增加了很多新功能&#xff0c;比如分区表性能优化…

skywalking9.4 链路追踪

下载&#xff0c;很慢很慢很慢&#xff01;&#xff01;&#xff01;&#xff01; jdk 使用jdk17 skywalking-apm 9.4 java-agent 9.0 idea 本地开发配置 第1行配置按实际来&#xff1b; 第2行自定义&#xff0c;一般和微服务名称相同&#xff1b; 第3行ip写安装的机器ip,端…

OpenGL3.3_C++_Windows(5)

变换 && 3D空间的2D图形 /\/\/\/\/\//\/\/\/\/\/\/\/\//\/\///\/\/\/\//\/\/\/\//\//\/\/\/\/\\/GLM库从0.9.9版本起&#xff0c;默认会将矩阵类型初始化为一个零矩阵&#xff08;所有元素均为0&#xff09;&#xff0c;而不是单位矩阵&#xff08;对角元素为1&#…

如何更新 iOS 18 Beta 版本?具体步骤总结

如何更新 iOS 18 Beta 想必有一些用户已经迫不及待的想要知道怎么更新 iOS 18 Beta 版本了吧&#xff0c;下面就给大家总结了具体的操作步骤&#xff1a; 在更新 iOS 18 Beta 版本之前记得我们需要将手机的数据进行备份&#xff0c;大家可以自行选用备份软件比如 iCloud 等。…

Eureka到Nacos迁移实战:解决配置冲突与启动异常

问题&#xff1a;Eureka到Nacos迁移实战&#xff1a;解决配置冲突与启动异常 在进行微服务架构升级&#xff0c;特别是注册中心从Eureka转向Nacos的过程中&#xff0c;我遇到了一个典型的技术挑战。目标是为了减少因配置变更导致的服务重启频率&#xff0c;我决定拥抱Nacos以其…

如何避免重复创建线程?创建线程池的方式有哪些?各自优缺点有哪些?

如何避免重复创建线程?创建线程池的方式有哪些?各自优缺点有哪些? 1、案例分析2、创建线程池的方式1、案例分析 我们在使用一些app的时候,应该都收到过消息推送,它们往往依赖消息推送服务实现。事实上,互联网大厂都有自己的消息推送服务(又名Message Push Server),通…

详解 Spring Security:全面保护 Java 应用程序的安全框架

详解 Spring Security&#xff1a;全面保护 Java 应用程序的安全框架 Spring Security 是一个功能强大且高度可定制的框架&#xff0c;用于保护基于 Java 的应用程序。它为身份验证、授权、防止跨站点请求伪造 (CSRF) 等安全需求提供了解决方案。下面将更详细地介绍 Spring Se…

复星杏脉算法面经2024年5月16日面试

复星杏脉算法面经2024年5月 面试记录&#xff1a;3个部分1. 自己介绍 2. 问八股 3.代码题先自我介绍20分钟问问题1. 梯度爆炸怎么解决&#xff0c;三个解决方案&#xff1a;梯度裁剪&#xff08;Gradient Clipping&#xff09;正则化&#xff08;Regularization&#xff09;调整…

Flink Sql:四种Join方式详解(基于flink1.15官方文档)

JOINs flink sql主要有四种连接方式&#xff0c;分别是Regular Joins、Interval Joins、Temporal Joins、lookup join 1、Regular Joins&#xff08;常规连接 &#xff09; 这种连接方式和hive sql中的join是一样的&#xff0c;包括inner join&#xff0c;left join&#xff…

java Springboot网上音乐商城(源码+sql+论文)

1.1 研究目的和意义 随着市场经济发展&#xff0c;尤其是我国加入WTO &#xff0c;融入经济全球化潮流&#xff0c;已进入国内外市场经济发展新时期&#xff0c;音乐与市场联系越来越紧密&#xff0c;我国音乐和网上业务也进入新历史发展阶段。为了更好地服务于市场&#xff0…

Studio One 6.6.2 for Mac怎么激活,有Studio One 6激活码吗?

如果您是一名音乐制作人&#xff0c;您是否曾经为了寻找一个合适的音频工作站而苦恼过&#xff1f;Studio One 6 for Mac是一款非常适合您的MacBook的音频工作站。它可以帮助您轻松地录制、编辑、混音和发布您的音乐作品。 Studio One 6.6.2 for Mac具有直观的界面和强大的功能…

全网爆火《pvz植物大战僵尸杂交版》最新安装包,Android、Windows、ios安装包+教程!

今天阿星想和大家分享一个最近在B站上引起轰动的老游戏——《植物大战僵尸》&#xff01; 是的&#xff0c;你没听错&#xff0c;就是那个曾经让我们熬夜到天亮&#xff0c;一关接一关挑战的游戏。 让我们来聊聊&#xff0c;这款游戏怎么就突然又火了起来呢&#xff1f; 原来…

AI玩具来了,它怎么样?

90后的我们&#xff0c;是AI时代的见证者。20后的小孩&#xff0c;才是AI时代的原著民。当ChatGPT们改变着大人的工作方式&#xff0c;我觉得&#xff0c;是时候让孩子们的玩具也更聪明些了吧。于是&#xff0c;在六一前夕&#xff0c;我用市面上的AI语音对话套件给娃DIY了一套…

简单的基于Transformer的滚动轴承故障诊断(Pytorch)

递归神经网络在很长一段时间内是序列转换任务的主导模型&#xff0c;其固有的序列本质阻碍了并行计算。因此&#xff0c;在2017年&#xff0c;谷歌的研究人员提出了一种新的用于序列转换任务的模型架构Transformer&#xff0c;它完全基于注意力机制建立输入与输出之间的全局依赖…