2012-12-23 69 views
11

Tôi đang cố gắng nắm bắt chức năng fft của Python, và một trong những điều kỳ lạ mà tôi đã vấp phải là Parseval's theorem dường như không áp dụng, vì nó mang lại sự khác biệt khoảng 50 bây giờ, trong khi nó phải được 0.Định lý của Parseval trong Python

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.fftpack as fftpack 

pi = np.pi 

tdata = np.arange(5999.)/300 
dt = tdata[1]-tdata[0] 

datay = np.sin(pi*tdata)+2*np.sin(pi*2*tdata) 
N = len(datay) 

fouriery = abs(fftpack.rfft(datay))/N 

freqs = fftpack.rfftfreq(len(datay), d=(tdata[1]-tdata[0])) 

df = freqs[1] - freqs[0] 

parceval = sum(datay**2)*dt - sum(fouriery**2)*df 
print parceval 

plt.plot(freqs, fouriery, 'b-') 
plt.xlim(0,3) 
plt.show() 

tôi khá chắc chắn rằng đó là một yếu tố bình thường, nhưng tôi dường như không thể tìm thấy nó, như tất cả các thông tin mà tôi có thể tìm thấy về chức năng này scipy.fftpack.rfft documentation.

Trả lời

15

Yếu tố chuẩn hóa của bạn đến từ cố gắng áp dụng định lý Parseval cho biến đổi Fourier của tín hiệu liên tục thành một chuỗi rời rạc. Trên bảng điều khiển bên của the wikipedia article on the Discrete Fourier transform có một số cuộc thảo luận về mối quan hệ của biến đổi Fourier, chuỗi Fourier, Biến đổi Fourier rời rạc và lấy mẫu với Dirac combs.

Để tạo một câu chuyện dài ngắn, Parseval's theorem, when applied to DFTs, không yêu cầu tích hợp, nhưng tổng kết: 2*pi bạn đang tạo bằng cách nhân số dtdf tóm tắt của bạn.

Cũng lưu ý rằng, vì bạn đang sử dụng scipy.fftpack.rfft, những gì bạn nhận được không đúng DFT của dữ liệu của bạn, nhưng chỉ là một nửa tích cực của nó, vì âm sẽ đối xứng với nó. Do bạn chỉ thêm một nửa dữ liệu, cộng với số 0 trong thuật ngữ DC, có thiếu 2 để truy cập vào số 4*pi mà @unutbu tìm thấy.

Trong mọi trường hợp, nếu datay giữ chuỗi, bạn có thể xác minh định lý Parseval như sau:

fouriery = fftpack.rfft(datay) 
N = len(datay) 
parseval_1 = np.sum(datay**2) 
parseval_2 = (fouriery[0]**2 + 2 * np.sum(fouriery[1:]**2))/N 
print parseval_1 - parseval_2 

Sử dụng scipy.fftpack.fft hoặc numpy.fft.fft tổng thứ hai không cần phải thực hiện trên một hình thức kỳ lạ như vậy:

fouriery_1 = fftpack.fft(datay) 
fouriery_2 = np.fft.fft(datay) 
N = len(datay) 
parseval_1 = np.sum(datay**2) 
parseval_2_1 = np.sum(np.abs(fouriery_1)**2)/N 
parseval_2_2 = np.sum(np.abs(fouriery_2)**2)/N 
print parseval_1 - parseval_2_1 
print parseval_1 - parseval_2_2 
+0

Cảm ơn bạn đã sửa! – unutbu

+0

Để lưu ý: một phần đây là khía cạnh của vấn đề chung mà các số dấu phẩy động không giống với số thực. – Marcin

+0

@Marcin Yep, đã thay đổi '==' thành '-' trong séc ... – Jaime