Nie zrobiłem DSP za jakiś czas, ale wręczamSpodziewałem się, że moje zrozumienie podstaw do zejścia tak daleko. Mam skrypt, w którym łączę ton ze złożonym wykładniczym. Rezultat, którego oczekuję, będzie przesunięty. Mój wynik jest dość nieoczekiwany - dostaję 3 tony i żaden z nich nie ma takiej częstotliwości, jakiej oczekuję. Czy ktoś może wyjaśnić, dlaczego otrzymuję te wyniki?
Oto skrypt.
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()
Odpowiedzi:
1 dla odpowiedzi № 1Korzystasz z właściwość przesunięcia częstotliwości. (Patrz np. http://ocw.usu.edu/Electrical_and_Computer_Engineering/Signals_and_Systems/5_6node6.html; przewiń w dół do sekcji oznaczonej Nieruchomość z przesunięciem częstotliwości.) To znaczy, jeśli transformata Fouriera f(t) jest F(w), następnie transformata Fouriera f(t)*exp(j*w0*t) jest F(w - w0). Ekspresja f(t)*exp(j*w0*t) jest punktowym mnożeniem f(t) i exp(j*w0*t), nie skręt.
Aby zobaczyć oczekiwany wynik, zastąp to:
res = scipy.signal.fftconvolve(f100, f200j, "full")
z
res = f100 * f200j
Wynik jest łatwiejszy do sprawdzenia, jeśli zmodyfikujesz funkcję wydruku w następujący sposób:
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)
i dodaj
from scipy.fftpack import fft, fftfreq
na górze skryptu.
Zobaczysz, że piki przy -8000 i 8000 na wykresie FFT f100 są przesunięte do -7000 i 9000 na wykresie FFT res.
0 dla odpowiedzi nr 2
To, czego chcesz, to mnożenie w dziedzinie częstotliwości (przesunięty ton), więc musisz to zrobić:
res = f100 * f200j
zamiast
res = scipy.signal.fftconvolve(f100, f200j, "full")
Poniżej znajduje się zdjęcie widma (nie FFT) input (czarny) i res (niebieski) sygnał. Przesunięty sygnał pokazuje dirak przy 9 kHz (8000 przesunięty w górę o 1000).
a poniżej jest moja funkcja 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)|")