In [305]:
import numpy as np
from math import *
import pandas as pd
from pylab import *

In [306]:
points = float(254)
modes = float(3)
seed = 5

In [307]:
def f(a, b, k1, k2, x, y):
    return a*cos(k1*x+k2*y) + b*sin(k1*x+k2*y)
def fx(a, b, k1, k2, x, y):
    return k1*(- a*sin(k1*x+k2*y) + b*cos(k1*x+k2&y) )
def fy(a, b, k1, k2, x, y):
    return k2*(- a*sin(k1*x+k2*y) + b*cos(k1*x+k2&y) )
def fxx(a, b, k1, k2, x, y):
    return k1*k1*(-a*cos(k1*x+k2*y) - b*sin(k1*x+k2&y))
def fyy(a, b, k1, k2, x, y):
    return k2*k2*(-a*cos(k1*x+k2*y) - b*sin(k1*x+k2&y))
def fxy(a, b, k1, k2, x, y):
    return k1*k2*(-a*cos(k1*x+k2*y) - b*sin(k1*x+k2&y))
def fxxx(a, b, k1, k2, x, y):
    return k1*k1*k1*(a*sin(k1*x+k2*y) - b*cos(k1*x+k2&y))
def fyyy(a, b, k1, k2, x, y):
    return k2*k2*k2*(a*sin(k1*x+k2*y) - b*cos(k1*x+k2&y))
def fxyy(a, b, k1, k2, x, y):
    return k1*k2*k2*(a*sin(k1*x+k2*y) - b*cos(k1*x+k2&y))
def fyxx(a, b, k1, k2, x, y):
    return k1*k1*k2*(a*sin(k1*x+k2*y) - b*cos(k1*x+k2&y))

Пример для трёх гармоник


In [315]:
mu, sigma = 0.0, 1.0 

a_coef = np.random.normal(mu, sigma, (2*int(points) + 1, int(points) + 1))
b_coef = np.random.normal(mu, sigma, (2*int(points) + 1, int(points) + 1))

In [316]:
def P(k1, k2, K0):
    if (k1 > K0) | (k2 > K0):
        return 0
    else:
        return 1

In [317]:
def func (i, j, a_coef, b_coef, modes, N):
    y = 0.0
    for k1 in range(0, 2*int(modes)): 
        for k2 in range(0, int(modes)): 
            y = y + (a_coef[k1][k2]*cos(2*pi*(i*(k1 - int(modes))+j*k2)/N) + b_coef[k1][k2]*sin(2*pi*(i*(k1 - int(modes))+j*k2)/N))
            #y = y + (a_coef[k1][k2]*cos(2*pi*(i*k1+j*k2)/N) + a_coef[k1][k2]*sin(2*pi*(i*k1+j*k2)/N))

    return y
#((float(k1 - int(modes)))**2-(float(k2))**2)**P(k1, k2, int(modes))*

In [318]:
file = open("f1.dat", 'w')
for j1 in xrange(0, int(points)):
    for j2 in xrange (0, int(points)):
        x = j1/(points)
        y = j2/(points)
        file.write('{}    {}    {}\n'.format(y, x, func(j1, j2, a_coef, b_coef, int(modes), points)))
file.close()

In [312]:
def func (i, j, a_coef, b_coef, modes, N):
    y = 0.0
    for k1 in range(0, 2*int(modes)): 
        for k2 in range(0, int(modes)): 
            y = y + (a_coef[k1][k2]*cos(2*pi*(i*(k1 - int(modes))+j*k2)/N) + b_coef[k1][k2]*sin(2*pi*(i*(k1 - int(modes))+j*k2)/N))
            #y = y + (a_coef[k1][k2]*cos(2*pi*(i*k1+j*k2)/N) + a_coef[k1][k2]*sin(2*pi*(i*k1+j*k2)/N))
    
    return y
#2*((float(k1 - int(modes)))*(float(k2)))*

In [319]:
mu, sigma = 0.0, 1.0 

a_coef = np.random.normal(mu, sigma, (2*int(points) + 1, int(points) + 1))
b_coef = np.random.normal(mu, sigma, (2*int(points) + 1, int(points) + 1))

In [320]:
file = open("f2.dat", 'w')
for j1 in xrange(0, int(points)):
    for j2 in xrange (0, int(points)):
        x = j1/(points)
        y = j2/(points)
        file.write('{}    {}    {}\n'.format(y, x, func(j1, j2, a_coef, b_coef, int(modes), points)))
file.close()

np.fft.fft2(x, 4)