| root2 = sqrt(2.0);N = 8,r = 2,sign = −1; |
| root_2N = sqrt(2.0 * N − 1.0);a = 2.0/root_2N; |
| for k = 0: N − 1; |
| if (mod(k, N)==0) |
| pk = 1.0/root2; |
| else pk = 1.0; |
| end; |
| M = N + N − 1; |
| theta = 2 * r * k * pi/M;Ar = cos(theta);Br = sin(theta); |
| sign = −sign; |
| p1 = 0.0;p2 = 0.0;p3 = 0.0;p4 = 0.0;p5 = 0.0 |
| for x = 0 : r − 1; |
| angle1 = (x + x + 1) * k * pi/M; |
| angle2 = (x + x) * k * pi/M; |
| sin 1 = sin(angle1); |
| sin 2 = sin(angle2); |
| cos 1 = cos(angle1); |
| cos 2 = cos(angle2); |
| if (x==0) |
| pk1 = 1.0/root2; |
| else |
| pk1 = 1.0; |
| end; |
| p1 = p1 + pk1 * (fnew(r − x)) * sin 2; |
| p2 = p2 + (fold(r − x)) * sin 1; |
| p3 = p3 + (−fold(2 * r − x) + Ar * fold(r − x)) * cos 1; |
| p4 = p4 + (fnew(2 * r − x) − pk1 * Ar * fnew(r − x)) * cos 2; |
| end; |
| p5 = (1.0/sqrt(2.0) − 1) * sign * (f(N) − fnew(r) * Ar + fnew(2 * r)); |
| term1 = a * pk * (−sign * Br * p1 + Br * p2 + p3 + sign * p4 + p5); |
| CNEW(k + 1) = 2 * Ar * C(k + 1) − Cold(k + 1) + term1 |
| end; |