有趣的傅里叶变换与小波变换对比(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;并取得了胜利。 迈…

课时157:脚本发布_简单脚本_功能函数

2.1.3 功能函数 学习目标 这一节&#xff0c;我们从 基础知识、简单实践、小结 三个方面来学习 基础知识 简介 需求&#xff1a;三条命令其实是一个组合&#xff0c;实现的是一个功能简单实践 实践 查看脚本内容 #!/bin/bash # 功能&#xff1a;打包代码 # 版本: v0.3 #…

代码随想录第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;我们将深入探讨监督学习的概念、机器学习的要求以及学习和提高预测准确…

Kotlin 协程:从基础概念到开发实践

前言 上一篇文章 深入理解Android多线程开发:场景应用与解决方案解析 针对Android开发中的多线程应用场景和相应的解决方案做了一个梳理。 总结出了Android开发中多线程编程的几个重要点: 资源复用和优化切线程任务编排并结合示例说明了Kotlin协程在处理上述问题时的优势。 …

如何删除:systemctl status ztncui

要删除 systemctl status ztncui&#xff0c;需要卸载与该服务相关的程序&#xff0c;并停止和禁用该服务。以下是详细的步骤&#xff1a; 停止服务&#xff1a; sudo systemctl stop ztncui 禁用服务&#xff1a; sudo systemctl disable ztncui 删除服务文件&#xff1a;…

归纳贪心好题

很有趣的一道归纳贪心题目 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…

图像分隔和深度成像技术为什么受市场欢迎-数字孪生技术和物联网智能汽车技术的大爆发?分析一下图像技术的前生后世

图像分隔和深度成像是计算机视觉和图像处理领域的两项重要技术&#xff0c;它们各自有不同的技术基础和要点。 图像分隔技术基础&#xff1a; 机器学习和模式识别&#xff1a; 图像分隔通常依赖于机器学习算法&#xff0c;如支持向量机&#xff08;SVM&#xff09;、随机森林…

算法题--华为od机试考试(开源项目热度榜单、API集群负载统计、分月饼)

目录 开源项目热度榜单 题目描述 输入描述 输出描述 示例1 输入 输出 解析 答案 API集群负载统计 题目描述 输入描述 输出描述 示例1 输入 输出 解析 答案 分月饼 题目描述 输入描述 输出描述 示例1 输入 输出 说明 示例2 输入 输出 说明 示例3…

Python 操作Redis

本篇介绍在Python中如何使用Redis。 目录 安装类库 引入类库 操作读写 使用pipeline 封装 封装redis 操作 调用 总结 安装类库 首先通过pip list命令查看&#xff0c; 本地安装类库中是否有redis类库&#xff0c; 如果没有开始安装。 安装命令如下&#xff1a; pip…

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

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

(简单)html转图片-html2canvas

官方文档 https://html2canvas.hertzen.com/ 案例代码 自己创建一个html然后试一下 <!DOCTYPE html> <html> <head><title>HTML to Image Example</title><script src"https://cdnjs.cloudflare.com/ajax/libs/html2canvas/1.3.2/ht…

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可以轻松实现这个过程。本文将介绍如…

Python 学习 第二册 第11章 文件

----用教授的方式学习 目录 11.1 打开文件 11.2 文件的基本方法 11.2.1 读取和写入 11.2.2 使用管道重定向输出 11.2.3 读取和写入行 11.2.4 关闭文件 11.3 迭代文件内容 11.3.1 每次一个字符&#xff08;或字节&#xff09; 11.3.2 每次一行 11.3.3 读取所有内容 …

make menuconfig | allyesconfig | allnoconfig【笔记】

make menuconfig, make allyesconfig, 和 make allnoconfig 是在编译Linux内核或某些其他使用类似配置系统的开源项目时使用的命令。这些命令用于生成内核配置文件&#xff0c;该文件决定了内核编译时哪些功能会被包含或排除。 下面是这三个命令的简要说明&#xff1a; 1、ma…

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;更别…