时域FFT

Question

我有下面的方波函数，只关注从-p0到p0的一个周期。我很难理解当我执行fft时到底发生了什么，我如何知道要使用哪些频率？

我手工计算了我的fft，应该类似于我在底部的函数。我做错了什么？

import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft,fftfreq,ifft
###################TIME/MOMENTUM DOMAIN########################
sample_rate = 1024
N = 2*sample_rate

p0 = 2
t = np.linspace(-p0,p0,N)
y = np.zeros_like(t)
y[(t>-p0)*(t<p0)] = np.sqrt(1/(2*p0)) 

#################Frequency/Position Domain#####################
frequency = np.linspace (0.0, 512, int (N/2)) #creates all the neccessary frequencies
fft_values = fft(y)
yf = 2/N * np.abs (fft_values [0:np.int (N/2)])
plt.plot(frequency,yf)
plt.show()




##############Calculated FFT########################
x = np.linspace(-p0,p0,1000)
y = [((np.sqrt(1/(np.pi*p0)))*np.sin(p0*t))/t for t in x]

Answer 1

你的问题我不太清楚。但看起来您是在尝试比较FFT的实践和理论方面。

让我们举一个简单的例子来理解这一点。我们有函数x1(t) = cos(2pi(0.5)t)*w(t)，窗口为10ms，余弦为0.5 cos。它看起来就像这样。

# Install this library if you needed, just do pip install pip install scikit-dsp-comm
# importing necessary libraries
%pylab inline
import sk_dsp_comm.sigsys as ss
import scipy.signal as signal
from IPython.display import Image, SVG

#Notebook configuration
pylab.rcParams['savefig.dpi'] = 100 # default 72
%config InlineBackend.figure_formats=['svg'] # SVG inline viewing

fs = 4 # sampling rate in kHz
t = arange(-5,5,1/fs)
tau = 1
f1 = 0.5; # Frequency component, in [kHz]
omega = 2*pi*f1;

x1 = cos(omega*t);
figure(figsize=(6,5))
subplot(211)
plot(t,x1.real,'b')
grid()
ylim([-1.1,1.1])
xlim([-5,5])
title(r'0.5KHz cosine over 10ms windows')
xlabel(r'Time (ms)')
ylabel(r'$x_0(t)$');

# FT Exact Plot
f,X1 = ss.ft_approx(x1,t,2000)
subplot(212)
plot(f,abs(X1),'b')
#plot(f,angle(X0))
grid()
xlim([-2,2])
title(r'Spectrum Magnitude of 0.5KHz cosine over 10ms windows')
xlabel(r'Frequency (kHz)')
ylabel(r'$|X_0e(f)|$');
tight_layout()

​

​

现在用0.5 this余弦的20ms窗口绘制相同的函数

t = arange(-10,10,1/fs)
x2 = cos(omega*t);
subplot(211)
plot(t,x2.real,'b')
grid()
ylim([-1.1,1.1])
xlim([-10,10])
title(r'0.5KHz cosine over 20ms windows')
xlabel(r'Time (ms)')
ylabel(r'$x_0(t)$');

# FT Exact Plot
f,X2 = ss.ft_approx(x2,t,2000)
subplot(212)
plot(f,abs(X2),'b')
#plot(f,angle(X0))
grid()
xlim([-2,2])
title(r'Spectrum Magnitude of 0.5KHz cosine over 20ms windows')
xlabel(r'Frequency (kHz)')
ylabel(r'$|X_0e(f)|$');
tight_layout()

​

​

现在，让我们介绍一下采样率和DFT中的2000个点。

fs = 4 # sampling rate in kHz
W = 5
t = arange(-5,5,1/fs)
x4 = W/pi*sinc(W/pi*t)
figure(figsize=(6,2))
plot(t,x4,'b')
grid()
# ylim([-1.1,1.1])
xlim([-5,5])
title(r'Time Domain: $x_4(t),\ W = 5$ Hz')
xlabel(r'Time (s)')
ylabel(r'$x_4(t)$');
f,X4 = ss.ft_approx(x4,t,2000)
figure(figsize=(6,2))
plot(f,abs(X4),'b')
grid()
title(r'Frequency Domain: $X_4(f)$')
xlim([-1,1])
xlabel(r'Frequency (Hz)')
ylabel(r'$|X_4(f)|$');
figure(figsize=(6,2))
plot(f,20*log10(abs(X4)),'b')
grid()
title(r'Frequency Domain: $X_4(f)$ in dB')
ylim([-50,5])
xlim([-1,1])
xlabel(r'Frequency (Hz)')
ylabel(r'$|X_4(f)|$ (dB)');

​

​

现在总结一下，这是一个简单的例子，其中单个正弦和它的FFT是这样表示的。

from scipy import fftpack
N = 10

t = np.linspace(0, 1, 500)
x = np.sin(49 * np.pi * t)

X = fftpack.fft(x)

f, (ax0, ax1) = plt.subplots(2, 1)

ax0.plot(x,'b');ax0.grid()
ax0.set_ylim(-1.1, 1.1)

ax1.plot(fftpack.fftfreq(len(t)), np.abs(X),'b');ax1.grid()
ax1.set_ylim(0, 190);

​

​

更多细节和Basic introduction可以在这里找到。