计算物理学，FFT分析

Question

我为计算任务解决了以下问题，但是我得到了非常糟糕的成绩（67%），我想了解如何正确地完成这些问题，特别是Q1.b和Q3。请尽可能详细，我真的很想理解我的错误。

生成数据（正弦函数）。使用FFT进行分析： a）三个具有恒定但不同频率的波的叠加 b）频率随时间变化的波 用适当的轴绘制图形、采样频率、振幅和功率谱。

使用Exercise 1a）中的三个波，但将它们更改为具有相同的频率、相位和振幅。随着随机高斯分布噪声逐渐增加，对每个波进行污染。 1）对三个噪声受污染的波的叠加执行FFT。分析并绘制输出。 2）使用高斯函数过滤信号，绘制“清洁”波形，并分析结果。得到的波形是否100%干净？请解释。

#1(b)

tmin = -2*pi
tmax - 2*pi
delta = 0.01
t = arange(tmin, tmax, delta)
y = sin(2.5*t*t)
plot(t, y, '-')
title('Figure 2: Plotting a wave whose frequency depends on time ')
xlabel('Time (s)')
ylabel('Y(t)')
show()

#b.2
Fs = 150.0;  # sampling rate
Ts = 1.0/Fs; # sampling interval
t = np.arange(0,1,Ts) # time vector

ff = 5;   # frequency of the signal
y = np.sin(2*np.pi*ff*t)

n = len(y) # length of the signal
k = np.arange(n)
T = n/Fs
frq = k/T # two sides frequency range
frq = frq[range(n/2)] # one side frequency range

Y = np.fft.fft(y)/n # fft computing and normalization
Y = Y[range(n/2)]

#Time vs. Amplitude
plot(t,y)
title('Figure 2: Time vs. Amplitude')
xlabel('Time')
ylabel('Amplitude')
plt.show()

#Amplitude Spectrum
plot(frq,abs(Y),'r')
title('Figure 2a: Amplitude Spectrum')
xlabel('Freq (Hz)')
ylabel('amplitude spectrum')
plt.show()


#Power Spectrum
plot(frq,abs(Y)**2,'r')
title('Figure 2b: Power Spectrum')
xlabel('Freq (Hz)')
ylabel('power spectrum')
plt.show()
#Exercise 3:

#part 1
t = np.linspace(-0.5*pi,0.5*pi,1000)

#contaminating our waves with successively increasing white noise
y_1 = sin(15*t) + np.random.normal(0,0.2*pi,1000)
y_2 = sin(15*t) + np.random.normal(0,0.3*pi,1000)
y_3 = sin(15*t) + np.random.normal(0,0.4*pi,1000)
y = y_1 + y_2 + y_3 # superposition of three contaminated waves


#Plotting the figure 
plot(t,y,'-')
title('A superposition of three waves contaminated with Gaussian Noise')
xlabel('Time (s)')
ylabel('Y(t)')
show()

delta = pi/1000.0
n = len(y)     ## calculate frequency in Hz
freq = fftfreq(n, delta)  # Computing the FFT


Freq = fftfreq(len(y), delta)  #Using Fast Fourier Transformation to         #calculate frequencies
N = len(Freq)
fr = Freq[1:len(Freq)/2.0] 
A = fft(y)
XF = A[1:len(A)/2.0]/float(len(A[1:len(A)/2.0]))


# Amplitude spectrum for contaminated waves
plt.plot(fr, abs(XF))   
title('Figure 3a : Amplitude spectrum with Gaussian Noise')
xlabel('frequency')
ylabel('Amplitude')
show()

# Power spectrum for contaminated waves
plt.plot(fr,abs(XF)**2)
title('Figure 3b: Power spectrum with Gaussian Noise')
xlabel('frequency(cycles/year)')
ylabel('Power')
show()

 # part 2
 F_v = exp(-(abs(freq)-2)**2/2*0.5**2)
 spectrum = A*F_v   #Applying the Gaussian Filter to clean our waves
 new_y = ifft(spectrum) #Computing the inverse FFT
 plot(t,new_y,'-')
 title('A superposition of three waves after Noise Filtering')
 xlabel('Time (s)')
 ylabel('Y(t)')
 show()

- ಠ_ಠ

---

欢迎来到 Stack Overflow。你问过批改者你哪里做错了吗？我们通常不会回答这样广泛的问题，比如“为什么我在这个复杂的多部分作业中得了低分？”投票关闭。 - Beta

最好的方法可能是礼貌地询问您的老师/助教问题出在哪里。 - pv.

我认为这个问题提得很好，而且（与任务偏离）的错误很容易看出来。我建议多次进行相同的任务，以真正理解FFT的概念，实函数的FFT会在正负频率上对称，这也是为什么我们只需要保留正频率。最重要的是要认识到频率间隔是时间范围的倒数，频率范围（负+正一起）是时间间隔的倒数。因此，采样定理在FFT提供的计算频率上完全得到满足。 - roadrunner66

Answer 1

以下是一些代码和图像，我们希望看到的就是这样的效果。为了展示所有的三个波和它们的总和，我在三个嘈杂波的总和的图例中偏离了情节。请注意，在嘈杂波的强度谱中，你看不到太多东西。对于这些情况，也可以绘制谱的对数（np.log），这样你就可以更好地看到噪音。
在最后一个图中，我绘制了高斯滤波器和频谱（不同大小）而没有重新缩放，只是为了展示滤波器的应用位置。它有效地是一个低通滤波器（通过较低的频率），通过与接近于零的数字相乘来消除更高频率的噪音。

import numpy as np
import matplotlib.pyplot as p
%matplotlib inline 

#1(b)
p.figure(figsize=(20,16))
p.subplot(431)
t = np.arange(0,10, 0.001)  #units in seconds
#cleaner to show the frequency change explicitly than y = sin(2.5*t*t)
f= 1+ t*0.1  # linear up chirp, i.e. frequency goes up  , frequency units in Hz (1/sec)
y = np.sin(2* np.pi* f* t)
p.plot(t, y, '-')
p.title('Figure 2: Plotting a wave whose frequency depends on time ')
p.xlabel('Time (s)')
p.ylabel('Y(t)')


#b.2
Fs = 150.0;  # sampling rate
Ts = 1.0/Fs; # sampling interval
t = np.arange(0,1,Ts) # time vector

ff = 5;   # frequency of the signal
y = np.sin(2*np.pi*ff*t)

n = len(y) # length of the signal
k = np.arange(n)   ## ok, the FFT has as many points in frequency space, as the original in time
T = n/Fs    ## correct ;  T=sampling time,   the total frequency range is 1/sample time
frq = k/T # two sided frequency range
frq = frq[range(n/2)] # one sided frequency range
Y = np.fft.fft(y)/n # fft computing and normalization
Y = Y[range(n/2)]

# Amplitude vs. Time
p.subplot(434)
p.plot(t,y)
p.title('y(t)') # Amplitude vs Time is commonly said, but strictly not true, the amplitude is unchanging
p.xlabel('Time')
p.ylabel('Amplitude')

#Amplitude Spectrum
p.subplot(435)
p.plot(frq,abs(Y),'r')
p.title('Figure 2a: Amplitude Spectrum')
p.xlabel('Freq (Hz)')
p.ylabel('amplitude spectrum')

#Power Spectrum
p.subplot(436)
p.plot(frq,abs(Y)**2,'r')
p.title('Figure 2b: Power Spectrum')
p.xlabel('Freq (Hz)')
p.ylabel('power spectrum')

#Exercise 3:

#part 1
t = np.linspace(-0.5*np.pi,0.5*np.pi,1000)

# #contaminating our waves with successively increasing white noise
y_1 = np.sin(15*t) + np.random.normal(0,0.1,1000)   # no need to get pi involved in this amplitude
y_2 = np.sin(15*t) + np.random.normal(0,0.2,1000)   
y_3 = np.sin(15*t) + np.random.normal(0,0.4,1000)   
y = y_1 + y_2 + y_3 # superposition of three contaminated waves


#Plotting the figure 
p.subplot(437)
p.plot(t,y_1+2,'-',lw=0.3)
p.plot(t,y_2,'-',lw=0.3)
p.plot(t,y_3-2,'-',lw=0.3)
p.plot(t,y-6 ,lw=1,color='black')
p.title('A superposition of three waves contaminated with  Gaussian Noise')
p.xlabel('Time (s)')
p.ylabel('Y(t)')


delta = np.pi/1000.0
n = len(y)     ## calculate frequency in Hz
# freq = np.fft(n, delta)  # Computing the FFT   <-- wrong, you don't calculate the FFT from a number, but from a time dep. vector/array
# Freq =  np.fftfreq(len(y), delta)  #Using Fast Fourier Transformation to         #calculate frequencies
# N = len(Freq)
# fr = Freq[1:len(Freq)/2.0] 
# A = fft(y)
# XF = A[1:len(A)/2.0]/float(len(A[1:len(A)/2.0]))   

# Why not do as before? 
k = np.arange(n)   ## ok, the FFT has as many points in frequency space, as the original in time
T = n/Fs    ## correct ;  T=sampling time,   the total frequency range is 1/sample time
frq = k/T # two sided frequency range
frq = frq[range(n/2)] # one sided frequency range
Y = np.fft.fft(y)/n # fft computing and normalization
Y = Y[range(n/2)]



# Amplitude spectrum for contaminated waves
p.subplot(438)
p.plot(frq, abs(Y))   
p.title('Figure 3a : Amplitude spectrum with Gaussian Noise')
p.xlabel('frequency')
p.ylabel('Amplitude')


# Power spectrum for contaminated waves
p.subplot(439)
p.plot(frq,abs(Y)**2)
p.title('Figure 3b: Power spectrum with Gaussian Noise')
p.xlabel('frequency(cycles/year)')
p.ylabel('Power')


# part 2

p.subplot(4,3,11)
F_v = np.exp(-(np.abs(frq)-2)**2/2*0.5**2) ## this is a Gaussian, plot it separately to see it; play with the values
cleaned_spectrum = Y*F_v   #Applying the Gaussian Filter to clean our waves ## multiplication in FreqDomain is convolution in time domain
p.plot(frq,F_v)
p.plot(frq,cleaned_spectrum)

p.subplot(4,3,10)  
new_y = np.fft.ifft(cleaned_spectrum) #Computing the inverse FFT of the cleaned spectrum to see the cleaned wave
p.plot(t[range(n/2)],new_y,'-')
p.title('A superposition of three waves after Noise Filtering')
p.xlabel('Time (s)')
p.ylabel('Y(t)')

- roadrunner66

---

你知道你有多厉害！！！非常感谢，这真的解决了很多问题，尤其是你在代码中的注释。 - ಠ_ಠ