Return to edu_fourier_0.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

Renaming:
 bench-1 --> bench-1.rr
 beta    --> beta.rr
 fourier-0 --> edu_fourier_0.rr
 hg21 --> taka_hg21.rr
 igraph --> igraph.rr
 sets --> oh_sets.rr
 uldet --> uldet.rr

/*  $OpenXM: OpenXM/src/asir-contrib/packages/src/edu_fourier_0.rr,v 1.1 2005/04/08 06:25:12 takayama Exp $ */
/* Old, fourier-0, see Attic */
/* 
   Name: edu_fourier_0
   Category: Education
   Subject: Signal Analysis. Digital Consine transformation.
   Main: asir.  Plugin: ox_sm1_gnuplot
   Documents: 
   See also: ~taka/this/Joho/fourier.c
*/
/*
   Sample Log:
      load("names.rr");
      gnuplot.start();
      load("./fourier-0");
      plotFourier_cos(2);  Displays graphs of the fourier transform of
                           cos(2*X) and the original cos(2*X).
      Change the value of NS. Then, you may get a better approximation.
      For example, compare NS=8 and NS=32.
      NS=8 is used in JPEG.
*/
if (ox_getenv("LANG") == "ja_JP.ujis") {
  print("$B;HMQNc(B 1:")$
  print("load(\"xm\")$")$
  print("gnuplot.start()$")$
  print("plotFourier_cos(2)$ --- cos(2*X) $B$N(B fourier $BE83+$N6a;w$H(B cos(2*X) $B$rI=<((B")$
/*  print("NS = 8 $B$G$"$k$,(B, $B$3$NCM$rJQ$($k$H(B fourier $BE83+$N78?t$N7W;;$N?tCM@QJ,$r$h$j$3$^$+$$(B mesh $B$G$*$3$J$&(B.")$ */
}else{
  print("Example 1:")$
  print("load(\"xm\")$")$
  print("gnuplot.start()$")$
  print("plotFourier_cos(2)$ --- cos(2*X) $B$N(B fourier $BE83+$N6a;w$H(B cos(2*X) $B$rI=<((B")$
}$


#define USE_GNU_PLOT 
/*
     F(x) = \sum_{k=1}^\infty a_k cos(k x).
     \frac{\pi}{2} a_k = \int_{0}^\pi F(x) cos(k x) dx
     a_0 = \int_{0}^\pi F(x) dx

     fourier(N,F) computes a_N for the rational N.
     fourier_cos_kx(N,K) computes a_N for cos(K*X).
*/

#define NS  8  /* Number of segment. [0,@pi] is divided into NS segments. */

/* N : number, F : Rational */
/*                You cannot use sin and cos in F. */
def fourier(N,F) {
  Sum = 0;
  for (I=0; I<NS; I++) {
    Sum = Sum+eval(subst(F,x,@pi*I/NS))*eval(cos(@pi*I*N/NS))*eval(@pi/NS);
  }
  if (N != 0) return(Sum*2/eval(@pi));
  else return(Sum/eval(@pi));
}

def invFourier(X,Table) {
  N = length(Table);
  Sum = 0;
  for (I=0; I<N; I++) {
    Sum = Sum+eval(Table[I]*cos(I*X));
  }
  return(Sum);
}


def fourier_cos_kx(N,K) {
  Sum = 0;
  for (I=0; I<NS; I++) {
    Sum = Sum+eval(cos(K*@pi*I/NS))*eval(cos(@pi*I*N/NS))*eval(@pi/NS);
  }
  if (N != 0) return(Sum*2/eval(@pi));
  else return(Sum);
}

/* Fourier expansion of cos(K*X) */
def plotFourier_cos(K) {
#ifdef USE_GNU_PLOT
  gnuplot.plot_dots([ ],0);
#endif
  G = [ ];
  Fcoeff = map(fourier_cos_kx,[0,1,2,3,4,5,6,7,8],K);
  print("Fourier coefficients : [a_0, a_1, ..., a_9] ------------------")$
  print(Fcoeff)$
  for (I=0; I<NS; I++) {
     X = eval(I*@pi/NS);
     G = append(G,[ [X,invFourier(X,Fcoeff)]]);
  }
  print("Values of \sum a_k cos(k*X) : [[x,value at x], ...] ----------")$
  print(G);
#ifdef USE_GNU_PLOT
  gnuplot.plot_dots(G,"lines");
  /*BUG: rtostr cannot convert bigfloat to a string.*/
#endif
  G2 = [ ];
  for (I=0; I<NS; I++) {
     X = eval(I*@pi/NS);
     G2 = append(G2,[ [X,eval(cos(K*X))]]);
  }
  print("Values of cos(P*X) : [[x,value at x], ...] ------------------")$
  print(G2);
#ifdef USE_GNU_PLOT
  gnuplot.plot_dots([ ],0);
  gnuplot.plot_dots(G2,"lines"); 
#endif
  G3 = [ ];
  for (I=0; I<NS; I++) {
     G3 = append(G3,[G[I][1]-G2[I][1]]);
  }
  print("The difference of y values.   -------------------------------")$
  print(G3);
}


end$