Làm thế nào để tính toán chuỗi Fourier trong Numpy?

Question

Tôi có một hàm tuần hoàn T và muốn biết cách lấy danh sách các hệ số Fourier. Tôi đã cố gắng sử dụng mô-đun fft từ numpy nhưng nó có vẻ chuyên dụng hơn để biến đổi Fourier hơn so với hàng loạt. Có thể nó thiếu kiến ​​thức toán học, nhưng tôi không thể thấy cách tính các hệ số Fourier từ fft.

Trợ giúp và/hoặc ví dụ được đánh giá cao.

Answer 1

Cuối cùng, điều đơn giản nhất (tính hệ số với một tổng Riemann) là/hiệu quả/cách mạnh mẽ cầm tay nhất để giải quyết vấn đề của tôi:

def cn(n): 
    c = y*np.exp(-1j*2*n*np.pi*time/period) 
    return c.sum()/c.size 

def f(x, Nh): 
    f = np.array([2*cn(i)*np.exp(1j*2*i*np.pi*x/period) for i in range(1,Nh+1)]) 
    return f.sum() 

y2 = np.array([f(t,50).real for t in time]) 

plot(time, y) 
plot(time, y2) 

mang lại cho tôi:

Answer 2

Bạn có danh sách các mẫu rời rạc về chức năng của mình hay chức năng của bạn có bị rời rạc không? Nếu vậy, Biến đổi Fourier rời rạc, được tính toán bằng thuật toán FFT, cung cấp các hệ số Fourier trực tiếp (see here).

Mặt khác, nếu bạn có biểu thức phân tích cho hàm, bạn có thể cần một trình giải toán toán học mang tính biểu tượng thuộc loại nào đó.

Answer 3

background information về thường trình là tài nguyên đầu tiên bạn nên tham khảo.

Answer 4

Numpy không phải là công cụ thích hợp thực sự để tính toán các thành phần chuỗi ô nhiễm bẩn, vì dữ liệu của bạn phải được lấy mẫu một cách rõ ràng. Bạn thực sự muốn sử dụng một cái gì đó như Mathematica hoặc nên sử dụng biến đổi fourier.

Để làm được điều đó, chúng ta hãy xem xét điều gì đó đơn giản là một sóng tam giác của giai đoạn 2pi, nơi chúng ta có thể dễ dàng tính toán hệ số Fourier (c_n = -i ((-1)^(n + 1))/n cho n > 0, ví dụ: c_n = {-i, i/2, -i/3, i/4, -i/5, i/6, ...} cho n = 1,2,3,4,5, 6 (sử dụng Sum (c_n exp (i 2 pi nx)) như Fourier series)

import numpy 
x = numpy.arange(0,2*numpy.pi, numpy.pi/1000) 
y = (x+numpy.pi/2) % numpy.pi - numpy.pi/2 
fourier_trans = numpy.fft.rfft(y)/1000 

Nếu bạn nhìn vào đầu tiên một số thành phần Fourier:.

array([ -3.14159265e-03 +0.00000000e+00j, 
     2.54994550e-16 -1.49956612e-16j, 
     3.14159265e-03 -9.99996710e-01j, 
     1.28143395e-16 +2.05163971e-16j, 
     -3.14159265e-03 +4.99993420e-01j, 
     5.28320925e-17 -2.74568926e-17j, 
     3.14159265e-03 -3.33323464e-01j, 
     7.73558750e-17 -3.41761974e-16j, 
     -3.14159265e-03 +2.49986840e-01j, 
     1.73758496e-16 +1.55882418e-17j, 
     3.14159265e-03 -1.99983550e-01j, 
     -1.74044469e-16 -1.22437710e-17j, 
     -3.14159265e-03 +1.66646927e-01j, 
     -1.02291982e-16 -2.05092972e-16j, 
     3.14159265e-03 -1.42834113e-01j, 
     1.96729377e-17 +5.35550532e-17j, 
     -3.14159265e-03 +1.24973680e-01j, 
     -7.50516717e-17 +3.33475329e-17j, 
     3.14159265e-03 -1.11081501e-01j, 
     -1.27900121e-16 -3.32193126e-17j, 
     -3.14159265e-03 +9.99670992e-02j, 

bỏ bê đầu tiên các thành phần mà gần 0 do độ chính xác của dấu phẩy động (~ 1e-16, bằng 0). phần khó khăn hơn là để thấy rằng các số 3.14159 (phát sinh trước khi chúng ta chia cho khoảng thời gian 1000) cũng phải được nhận biết là không, vì hàm này là định kỳ). Vì vậy, nếu chúng ta bỏ qua những hai yếu tố chúng tôi nhận được:

fourier_trans = [0,0,-i,0,i/2,0,-i/3,0,i/4,0,-i/5,0,-i/6, ... 

và bạn sẽ nhìn thấy số loạt Fourier đưa ra như mọi số khác (tôi đã không điều tra; nhưng tôi tin rằng các thành phần tương ứng với [c0, c- 1, c1, c-2, c2, ...]). Tôi đang sử dụng các quy ước theo wiki: http://en.wikipedia.org/wiki/Fourier_series.

Một lần nữa, tôi khuyên bạn nên sử dụng toán học hoặc hệ thống đại số máy tính có khả năng tích hợp và xử lý các chức năng liên tục.

Answer 5

Như các câu trả lời khác đã đề cập, có vẻ như những gì bạn đang tìm kiếm là một gói tính toán biểu tượng, do đó, không có vón cục là không phù hợp. Nếu bạn muốn sử dụng giải pháp dựa trên python miễn phí, thì sympy hoặc sage phải đáp ứng nhu cầu của bạn.

Answer 6

Đây là một câu hỏi cũ, nhưng vì tôi phải viết mã này, tôi đăng ở đây giải pháp sử dụng mô-đun numpy.fft, có khả năng nhanh hơn các giải pháp thủ công khác.

DFT là công cụ thích hợp để tính toán chính xác số các hệ số của chuỗi Fourier của hàm, được định nghĩa là biểu thức phân tích đối số hoặc hàm nội suy số trên một số điểm rời rạc.

này là việc thực hiện, cho phép để tính toán các hệ số thực có giá trị của chuỗi Fourier, hoặc các hệ số phức tạp có giá trị, bằng cách thông qua một thích hợp return_complex:

def fourier_series_coeff_numpy(f, T, N, return_complex=False): 
    """Calculates the first 2*N+1 Fourier series coeff. of a periodic function. 

    Given a periodic, function f(t) with period T, this function returns the 
    coefficients a0, {a1,a2,...},{b1,b2,...} such that: 

    f(t) ~= a0/2+ sum_{k=1}^{N} (a_k*cos(2*pi*k*t/T) + b_k*sin(2*pi*k*t/T)) 

    If return_complex is set to True, it returns instead the coefficients 
    {c0,c1,c2,...} 
    such that: 

    f(t) ~= sum_{k=-N}^{N} c_k * exp(i*2*pi*k*t/T) 

    where we define c_{-n} = complex_conjugate(c_{n}) 

    Refer to wikipedia for the relation between the real-valued and complex 
    valued coeffs at http://en.wikipedia.org/wiki/Fourier_series. 

    Parameters 
    ---------- 
    f : the periodic function, a callable like f(t) 
    T : the period of the function f, so that f(0)==f(T) 
    N_max : the function will return the first N_max + 1 Fourier coeff. 

    Returns 
    ------- 
    if return_complex == False, the function returns: 

    a0 : float 
    a,b : numpy float arrays describing respectively the cosine and sine coeff. 

    if return_complex == True, the function returns: 

    c : numpy 1-dimensional complex-valued array of size N+1 

    """ 
    # From Shanon theoreom we must use a sampling freq. larger than the maximum 
    # frequency you want to catch in the signal. 
    f_sample = 2 * N 
    # we also need to use an integer sampling frequency, or the 
    # points will not be equispaced between 0 and 1. We then add +2 to f_sample 
    t, dt = np.linspace(0, T, f_sample + 2, endpoint=False, retstep=True) 

    y = np.fft.rfft(f(t))/t.size 

    if return_complex: 
     return y 
    else: 
     y *= 2 
     return y[0].real, y[1:-1].real, -y[1:-1].imag 

Đây là một ví dụ về việc sử dụng:

from numpy import ones_like, cos, pi, sin, allclose 
T = 1.5 # any real number 

def f(t): 
    """example of periodic function in [0,T]""" 
    n1, n2, n3 = 1., 4., 7. # in Hz, or nondimensional for the matter. 
    a0, a1, b4, a7 = 4., 2., -1., -3 
    return a0/2 * ones_like(t) + a1 * cos(2 * pi * n1 * t/T) + b4 * sin(
     2 * pi * n2 * t/T) + a7 * cos(2 * pi * n3 * t/T) 


N_chosen = 10 
a0, a, b = fourier_series_coeff_numpy(f, T, N_chosen) 

# we have as expected that 
assert allclose(a0, 4) 
assert allclose(a, [2, 0, 0, 0, 0, 0, -3, 0, 0, 0]) 
assert allclose(b, [0, 0, 0, -1, 0, 0, 0, 0, 0, 0]) 

Và những âm mưu của các kết quả a0,a1,...,a10,b1,b2,...,b10 hệ số:

Đây là một thử nghiệm tùy chọn cho hàm, cho cả hai chế độ hoạt động. Bạn nên chạy sau ví dụ này hoặc xác định hàm định kỳ f và một khoảng thời gian T trước khi chạy mã.

# #### test that it works with real coefficients: 
from numpy import linspace, allclose, cos, sin, ones_like, exp, pi, \ 
    complex64, zeros 


def series_real_coeff(a0, a, b, t, T): 
    """calculates the Fourier series with period T at times t, 
     from the real coeff. a0,a,b""" 
    tmp = ones_like(t) * a0/2. 
    for k, (ak, bk) in enumerate(zip(a, b)): 
     tmp += ak * cos(2 * pi * (k + 1) * t/T) + bk * sin(
      2 * pi * (k + 1) * t/T) 
    return tmp 


t = linspace(0, T, 100) 
f_values = f(t) 
a0, a, b = fourier_series_coeff_numpy(f, T, 52) 
# construct the series: 
f_series_values = series_real_coeff(a0, a, b, t, T) 
# check that the series and the original function match to numerical precision: 
assert allclose(f_series_values, f_values, atol=1e-6) 

# #### test similarly that it works with complex coefficients: 

def series_complex_coeff(c, t, T): 
    """calculates the Fourier series with period T at times t, 
     from the complex coeff. c""" 
    tmp = zeros((t.size), dtype=complex64) 
    for k, ck in enumerate(c): 
     # sum from 0 to +N 
     tmp += ck * exp(2j * pi * k * t/T) 
     # sum from -N to -1 
     if k != 0: 
      tmp += ck.conjugate() * exp(-2j * pi * k * t/T) 
    return tmp.real 

f_values = f(t) 
c = fourier_series_coeff_numpy(f, T, 7, return_complex=True) 
f_series_values = series_complex_coeff(c, t, T) 
assert allclose(f_series_values, f_values, atol=1e-6)