从numpy简单卷积产生的意外结果

我在一段时间内还没有完成DSP，但是我并没有预料到我对基础知识的理解会变得如此糟糕。我有一个脚本，我正在用一个复杂的指数函数来调和一个音调。我期望的结果是一个移调。我的结果是非常意外的 - 我得到3声调，没有一个达到我期望的频率。有人能解释为什么我会得到这些结果吗？从numpy简单卷积产生的意外结果

这是脚本。

import sys 
import numpy 
import math 
import scipy 
from pylab import * 

def gen_tone(f, fs, length): 
    t = linspace(0, length, length * fs) 
    return cos(2.0 * pi * f * t) 

def gen_exp(f, fs, length): 
    t = linspace(0, length, length * fs) 
    return numpy.exp(1.0j * 2 * pi * f * t) 

def plot_fft(f, fs): 
    FFT = abs(scipy.fft(f, 1024))/f.size 
    figure() 
    plot(FFT) 

f100 = gen_tone(8000, 44100, 1) 
f200j = gen_exp(1000, 44100, 1) 

res = scipy.signal.fftconvolve(f100, f200j, 'full') 

plot_fft(f100, 44100) 
plot_fft(f200j, 44100) 
plot_fft(res, 44100) 

show()

来源

2014-02-07 Daniel

你期望结果如何（我的意思是说）？你有什么作为频率？ – usethedeathstar

正如沃伦的回答所说：你卷积了你的函数，而不是乘以它们。乘以复指数 - >频移。 – David

是的 - 这是我的失败。 – Daniel

您正在使用频移特性。（参见，例如http://ocw.usu.edu/Electrical_and_Computer_Engineering/Signals_and_Systems/5_6node6.html;稍微向下滚动到标记为移频属性部分）即，如果傅里叶变换的f(t)是F(w)，那么傅立叶变换的f(t)*exp(j*w0*t)是F(w - w0)。表达式f(t)*exp(j*w0*t)是f(t)和exp(j*w0*t),而不是卷积的逐点乘法。

一看就知道你所期望的结果，替换此：

res = scipy.signal.fftconvolve(f100, f200j, 'full')

与

res = f100 * f200j

结果是更容易地看到，如果你修改你的绘图功能如下：

def plot_fft(f, fs): 
    FFT = abs(fft(f, 1024))/f.size 
    freq = fftfreq(1024, 1.0/fs) 
    ndx = freq.argsort() 
    figure() 
    plot(freq[ndx], FFT[ndx]) 
    grid(True)

并添加

from scipy.fftpack import fft, fftfreq

位于脚本的顶部。

您将会看到f100的FFT图中-8000和8000处的峰在res的FFT图中转换为-7000和9000。

来源

2014-02-07 19:44:19

那么我做的操作等同于时间转换信号？我正在考虑时间平移规则。它在频域中乘以等于时域卷积的复指数。 – Daniel

时域中的卷积对应于频域中的逐点乘法。这就是为什么你有三个峰值，分别是-8000,8000和1000.如果你在上面提到的链接中向上滚动一下，你会看到“时移属性”的描述，它表示点对点乘法由复指数exp（-j * w * t0）给出的频域在时域中给出了't0'的偏移。 –

太棒了。感谢您的解释。 – Daniel

你想要的是在频域（移位音）的乘积，所以你必须这样做：

res = f100 * f200j

代替

res = scipy.signal.fftconvolve(f100, f200j, 'full')

下面的图片频谱（而不是FFT）的input（黑色）和res（蓝色）信号。移位的信号在9kHz处显示dirac（8000移位1000）。 enter image description here

这里下面是我plotSpectrum功能

def plotSpectrum(y,Fs): 
    """ 
    Plots a Single-Sided Amplitude Spectrum of y(t) 
    """ 
    n = len(y) # length of the signal 
    k = arange(n) 
    T = n/Fs 
    frq = k/T # two sides frequency range 
    frq = frq[arange(int(n/2))] # one side frequency range 

    Y = fft(y)/n # fft computing and normalization 
    Y = Y[arange(int(n/2))] 

    plot(frq,abs(Y),'r') # plotting the spectrum 
    xlabel('Freq (Hz)') 
    ylabel('|Y(freq)|')

来源

2014-02-07 13:31:11 RawBean

从numpy简单卷积产生的意外结果

回答

相关问题