從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簡單卷積產生的意外結果

回答

相關問題