使用numpy.fft进行卷积导致时间偏移

问题描述 投票:0回答:1

我想写一个完全在谱域的卷积代码。我在时间上取一个n个样本的尖峰序列(为了简单起见,下面的例子只有一个尖峰),然后用以下方法计算傅里叶序列。numpy.fft.fft. 我创建了一个m个样本的 "Ricker小波"(m << n),并使用以下方法计算其傅里叶序列 numpy.fft.fft但规定其输出的傅里叶序列为n个样本长。尖峰序列和小波都有相同的采样间隔。由此产生的卷积序列是偏移的(小波的峰值相对于尖峰沿时间轴偏移),这种偏移似乎取决于小波的大小,m。这种偏移似乎取决于小波的大小,m。

我认为这与以下参数有关。numpy.fft.fft(a, n=None, axis=-1, norm=None)尤其是 "轴 "参数。但是,我完全不理解这个参数的文档。

有谁能帮我理解为什么我会出现这种偏移(如果不清楚,让我明确地说,卷积序列中的小波峰值必须是输入尖峰序列中尖峰的同时采样)?

我的代码如下。

################################################################################
#
# import libraries
#
import math
import numpy as np
import scipy
import matplotlib.pyplot as plt
import os
from matplotlib.ticker import MultipleLocator
from random import random
# Define lists
#
Time=[]; Ricker=[]; freq=25; rickersize=51; timeiter=0.002; serieslength=501; TIMElong=[]; Reflectivity=[];
Series=[]; IMPEDANCE=[]; CONVOLUTION=[];
#
# Create ricker wavelet and its time sequence
#
for i in range(0,rickersize):
    time=(float(i-rickersize//2)*timeiter)
    ricker=(1-2*math.pi*math.pi*freq*freq*time*time)*math.exp(-1*math.pi*math.pi*freq*freq*time*time)
    Time.append(time)
    Ricker.append(ricker)
#
# Do various FFT operations on the Ricker wavelet:
#   Normal FFT, FFT of longer Ricker, Amplitude of the FFTs, their inverse FFTs and their frequency sequence
#
FFT=np.fft.fft(Ricker); FFTlong=np.fft.fft(Ricker,n=serieslength,axis=0,norm=None);
AMP=abs(FFT); AMPlong=abs(FFTlong);
RICKER=np.fft.ifft(FFT); RICKERlong=np.fft.ifft(FFTlong);
FREQ=np.fft.fftfreq(len(Ricker),d=timeiter); FREQlong=np.fft.fftfreq(len(RICKERlong),d=timeiter)
PHASE=np.angle(FFT); PHASElong=np.angle(FFTlong);
#
# Create a single spike in the otherwise empty (0) series of length 'serieslength' (=len(RICKERlong)
#   this spikes mimics a very simple seismic reflectivity series in time
#
for i in range(0,serieslength):
    time=(float(i)*timeiter)
    TIMElong.append(time)
    if i==int(serieslength/2):
        Series.append(1)
    else:
        Series.append(0)
#
# Do various FFT operations on the spike series
#   Normal FFT, Amplitude of the FFT, its inverse FFT and frequency sequence
#
FFTSeries=np.fft.fft(Series)
AMPSeries=abs(FFTSeries)
SERIES=np.fft.ifft(FFTSeries)
FREQSeries=np.fft.fftfreq(len(Series),d=timeiter)
#
# Do convolution of the spike series with the (long) Ricker wavelet in the frequency domain and see result via inverse FFT
#
FFTConvolution=[FFTlong[i]*FFTSeries[i] for i in range(len(Series))]
CON=np.fft.ifft(FFTConvolution)
CONVOLUTION=[CON[i].real for i in range(len(Series))]
#
# plotting routines
#
fig,axs = plt.subplots(nrows=1,ncols=3, figsize=(14,8))
axs[0].barh(TIMElong,Series,height=0.005, color='black')
axs[1].plot(Ricker,Time,color='black', linestyle='solid',linewidth=1)
axs[2].plot(CONVOLUTION,TIMElong,color='black', linestyle='solid',linewidth=1)
#
axs[0].set_aspect(aspect=8); axs[0].set_title('Reflectivity',fontsize=12); axs[0].yaxis.grid(); axs[0].xaxis.grid(); 
axs[0].set_xlim(-2,2); axs[0].set_ylim(min(TIMElong),max(TIMElong)); axs[0].invert_yaxis(); axs[0].tick_params(axis='both',which='major',labelsize=12);
#
axs[1].set_aspect(aspect=6.2); axs[1].set_title('Ricker',fontsize=12); axs[1].yaxis.grid(); axs[1].xaxis.grid(); 
axs[1].set_xlim(-1.0,1.02); axs[1].set_ylim(min(Time),max(Time)); axs[1].invert_yaxis(); axs[1].tick_params(axis='both',which='major',labelsize=12);
#
axs[2].set_aspect(aspect=8); axs[2].set_title('Convolution',fontsize=12); axs[2].yaxis.grid(); axs[2].xaxis.grid();
axs[2].set_xlim(-2,2); axs[2].set_ylim(min(TIMElong),max(TIMElong)); axs[2].invert_yaxis(); axs[2].tick_params(axis='both',which='major',labelsize=12);
#
fig.tight_layout()
fig.show()
####
python numpy fft convolution
1个回答
0
投票

结果据我所知,我的问题与python和numpy的特殊性无关。问题是 "循环卷积"。也就是说,两个数据序列的卷积是由两个序列的长度组合而成的,比较长。这必须在fft和ifft中加以说明。我不是这样做的。我仍然不知道具体如何处理这个问题,但现在我知道问题所在了,应该更简单了。

向那些试图回答我的畸形问题的人道歉。

© www.soinside.com 2019 - 2024. All rights reserved.