有趣的傅里叶变换与小波变换对比(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.mzph.cn/pingmian/28242.shtml

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

相关文章

摄影师在人工智能竞赛中与机器较量并获胜

摄影师在人工智能竞赛中与机器较量并获胜 自从生成式人工智能出现以来&#xff0c;由来已久的人机大战显然呈现出一边倒的态势。但是有一位摄影师&#xff0c;一心想证明用人眼拍摄的照片是有道理的&#xff0c;他向算法驱动的竞争对手发起了挑战&#xff0c;并取得了胜利。 迈…

代码随想录第29天|贪心算法part3

134.加油站 首先如果总油量减去总消耗大于等于零那么一定可以跑完一圈 每个加油站的剩余量rest[i]为gas[i] - cost[i] 从0开始累加rest[i]&#xff0c;和记为curSum&#xff0c;一旦curSum小于零&#xff0c;说明[0, i]区间都不能作为起始位置 因为我们一直维护的是一个剩余量大…

HCIA11 网络安全之本地 AAA 配置实验

AAA 提供 Authentication&#xff08;认证&#xff09;、Authorization&#xff08;授权&#xff09;和 Accounting&#xff08;计费&#xff09;三种安全功能。 • 认证&#xff1a;验证用户是否可以获得网络访问权。 • 授权&#xff1a;授权用户可以使用哪些服务。 •…

机器学习中的监督学习介绍

In this post well go into the concept of supervised learning, the requirements for machines to learn, and the process of learning and enhancing prediction accuracy. 在这篇文章中&#xff0c;我们将深入探讨监督学习的概念、机器学习的要求以及学习和提高预测准确…

归纳贪心好题

很有趣的一道归纳贪心题目 class Solution { public:int minimumAddedCoins(vector<int>& coins, int target) {sort(coins.begin(),coins.end());int n coins.size();int s 0,i0;int res 0;while(s<target){if(i<n&&coins[i]<s1)scoins[i];els…

顶顶通呼叫中心中间件-限制最大通话时间(mod_cti基于FreeSWITCH)

顶顶通呼叫中心中间件-限制最大通话时间(mod_cti基于FreeSWITCH) 一、最大通话时间 1、配置拨号方案 1、点击拨号方案 ->2、在框中输入通话最大时长->3、点击添加->4、根据图中配置->5、勾选continue。修改拨号方案需要等待一分钟即可生效 action"sched…

9M高速USB转接芯片CH347转双串口转I2C转SPI转JTAG转SWD

1、概述 CH347 TSSOP20封装和丝印 CH347 是一款高速 USB 总线转接芯片&#xff0c;通过 USB 总线提供异步串口、I2C 同步串行接口、SPI 同步串行接口和 JTAG 接口等。 在异步串口方式下&#xff0c;CH347 提供了 2 个高速串口&#xff0c;支持 RS485 串口收发使能控制、硬件流控…

LeetCode | 387.字符串中的第一个唯一字符

这道题可以用字典解决&#xff0c;只需要2次遍历字符串&#xff0c;第一次遍历字符串&#xff0c;记录每个字符出现的次数&#xff0c;第二次返回第一个出现次数为1的字符的下标&#xff0c;若找不到则返回-1 class Solution(object):def firstUniqChar(self, s):""…

Python自动化办公(一) —— 根据PDF文件批量创建Word文档

Python自动化办公&#xff08;一&#xff09; —— 根据PDF文件批量创建Word文档 在日常办公中&#xff0c;我们经常需要根据现有的PDF文件批量创建Word文档。手动操作不仅费时费力&#xff0c;而且容易出错。幸运的是&#xff0c;使用Python可以轻松实现这个过程。本文将介绍如…

linux的repo工具的入门

repo 是一个工具&#xff0c;用于管理 Git 仓库的集合&#xff0c;尤其在 Android 开发中被广泛使用。它是 Google 为 Android 项目开发的&#xff0c;以简化对大量 Git 仓库的管理。 主要特点 多仓库管理&#xff1a;repo 允许同时管理多个 Git 仓库&#xff0c;可以轻松执行…

MyBatis 的多级缓存机制是怎么样运作的?

引言&#xff1a;上周三&#xff0c;小 X 去面试一家中厂&#xff0c;其中面试官问到 MyBatis 的多级缓存机制是怎么样运行的&#xff1f;这个问题可以好好准备一下&#xff0c;很多人可能只会用 MyBatisPlus&#xff0c;简单的多表联查 SQL 语句可能都写不出来&#xff0c;更别…

Python 项目应该放弃requirements.txt?揭秘PDM的强大功能

目录 requirements.txt的局限性 PDM 的优势 如何使用 PDM 安装 PDM 初始化项目 添加依赖 管理依赖 示例代码 初始化项目并添加依赖 编写简单的 Flask 应用 运行应用 PDM高级功能 多环境管理 脚本管理 发布包 在 Python 项目中管理依赖项&#xff0c;最常见的方式…

Android APP memory统计方法

目录 进程的内存信息概述 关键的术语 测试步骤 测试步骤 数据处理 数据分析&#xff1a; 进程内存信息 Dumpsys meminfo -a PID Procrank Procmem PID 特殊内存信息 Mali ION(multi-media&#xff0c;gralloc) 进程地址空间信息 /proc/pid/smaps Showmap PID …

随笔-来了,安了

依照领导定的规矩&#xff0c;周五又去了分公司&#xff0c;赋能一线去了。到了地方就是开会->现场解决问题->干饭->开会过需求、提供解决方案&#xff0c;充实得厉害。强度也不小&#xff0c;中午干的一大碗饭&#xff0c;到五点就饿了。 六点带着分公司催着上线的需…

5000天后的世界

为何可以预见未来 1993年&#xff0c;在互联网的黎明时代&#xff0c;凯文凯利创办了《连线》杂志。他曾经采访过以比尔盖茨、史蒂夫乔布斯、杰夫贝佐斯为代表的一众风云创业家。《连线》杂志是全球发行的世界著名杂志&#xff0c;一直致力于报道科学技术带来的经济、社会变革…

【0基础学爬虫】爬虫基础之自动化工具 DrissionPage 的使用

概述 前三期文章中已经介绍到了 Selenium 与 Playwright 、Pyppeteer 的使用方法&#xff0c;它们的功能都非常强大。而本期要讲的 DrissionPage 更为独特&#xff0c;强大&#xff0c;而且使用更为方便&#xff0c;目前检测少&#xff0c;强烈推荐&#xff01;&#xff01;&a…

Google Earth Engine(GEE)——计算闪闪红星的ndvi的值和折线图(时序分析)

函数: ui.Chart.image.doySeries(imageCollection, region, regionReducer, scale, yearReducer, startDay, endDay)

手把手教你改造Sentinel Dashboard 实现配置持久化

一. 概述 Sentinel客户端默认情况下接收到 Dashboard 推送的规则配置后&#xff0c;可以实时生效。但是有一个致命缺陷&#xff0c;Dashboard和业务服务并没有持久化这些配置&#xff0c;当业务服务重启后&#xff0c;这些规则配置将全部丢失。 Sentinel 提供两种方式修改规则…

政务云参考技术架构

行业优势 总体架构 政务云平台技术框架图&#xff0c;由机房环境、基础设施层、支撑软件层及业务应用层组成&#xff0c;在运维、安全和运营体系的保障下&#xff0c;为政务云使用单位提供统一服务支撑。 功能架构 标准双区隔离 参照国家电子政务规范&#xff0c;打造符合标准的…

B3981 [信息与未来 2024] 图灵完备

题目描述 &#xff08;你不需要看懂这张图片&#xff1b;但如果你看懂了&#xff0c;会觉得它很有趣。&#xff09; JavaScript 是一种功能强大且灵活的编程语言&#xff0c;也是现代 Web 开发的三大支柱之一 (另外两个是 HTML 和 CSS)。灵活的 JavaScript 包含“自动类型转换…