C言語での矩形近似によるフーリエ級数展開が上手くいきません...

Question

前提・実現したいこと

次のような三角波、正弦波の矩形近似によるフーリエ級数展開を考えています。コメント文で示した部分(3か所)をどう記述すればよいかご教示いただければ幸いです。

f(t) = 1/2(1 - 2t/T) -2/T < t < T/2
1/2 t = ± T/2
f(t) = 1/2{ 1 + cos(2πt / T) -T/2 ≦ t ≦ T/2

該当のソースコード

#include <stdio.h>
#include <math.h>
#define M_PI 3.14159265 
#define MM 300   //フーリエ級数展開の打ち切り次数
#define Tp 1.0   //周期波形の周期

double aa[MM+1], bb[MM+1] ; //フーリエ展開係数 
int NN = 3 * MM ; //数値積分の分割数

double ft(double t)
{ 
  double zero = 1.0e-8 ;
  while ( t < -Tp/2.0 ){ t = t+Tp ;if( t > -Tp/2.0 ) break ;}
  while ( t >  Tp/2.0 ){ t = t-Tp ;if( t <  Tp/2.0 ) break ;}


   if( fabs( t+Tp/2.0 ) < zero || fabs( t-Tp/2.0) < zero ) return 0.5 ;
   else return 0.5 - t/Tp; //ここは間違っている可能性あり①


   return 0.5 + 0.5 * cos((2.0 * M_PI * t) / Tp); //ここは間違っている可能性あり②
}

double ak_dt(int k, double t)
{ 
  return ft(t) * cos( (double)k * 2.0 * M_PI / Tp * t ) ;
}

double fcoeff_ak( int k )
{ /*　関数f(t)のフーリエ展開級数 a_kを数値積分で求める */
        double t, i, ak = 0.0, dt = Tp/(double)NN ; 
        for( i = 0 ; i < NN ; i ++ ) {
          t = -Tp/2.0 + (double)i * dt ;
          ak = ak + sin(t) /  t ;  // ここがよくわかっていません。③
        }

        return ak ;
}

double bk_dt(int k, double t)
{ 
  return ft(t) * sin( (double)k * 2.0 * M_PI / Tp * t ) ;
}

double fcoeff_bk(int k)
{ 
        double t, i, bk = 0.0, dt = Tp/(double)NN ; 
        for( i = 0 ; i < NN ; i ++ ) {
          t = -Tp/2.0 + (double)i * dt ;
          bk = bk + sin(t) / t ;  // 同様にここがわかりません。④
        }
        return bk ;
}

double fxfreq( double *aa, double *bb, int M, double t)
{ 
        int k ; double ff = 0.0  ;
        for ( k = 0 ; k <= MM ; k++ ) ff = *aa*cos(((double)k * 2.0 * M_PI * t)/ Tp ) + *bb*sin(((double)k * 2.0 * M_PI * t)/Tp ) ;//ここは合っているのかどうかわかりません。⑤
        return ff;
}

int main()
{
  int i, k, n=2 ; double t, fh ;
  int KK = 1000 ;
  
  aa[0] = fcoeff_ak(0)/2.0 ; bb[0] = 0.0 ; 
  for (k=1; k<=MM; k++) {
    aa[k] = fcoeff_ak(k) ; bb[k] = fcoeff_bk(k) ; 
  }
  
  
  for (i=0; i<= n*KK; i++) {
    t = -Tp*(double)n/2.0 + Tp/(double)KK * (double)i ;
    fh = fxfreq( aa, bb, MM, t ) ;
    printf(" %10.8f   %10.8f   %10.8f\n", t, ft(t), fh );
  }
  return 0;
}

試したこと

実行してみましたが、正しい結果が得られませんでした。実行結果の一部分を掲載いたします。ft(t)は正しいと思います。fhもft(t)と近い値になるはずだと思うのですが、同じ値を行ったり来たりしています。
-0.00900000 0.50900000 -137.14071854
-0.00800000 0.50800000 -359.03903645
-0.00700000 0.50700000 359.03904488
-0.00600000 0.50600000 137.14070490
-0.00500000 0.50500000 -443.79666119
-0.00400000 0.50400000 137.14071399
-0.00300000 0.50300000 359.03903926
-0.00200000 0.50200000 -359.03904207
-0.00100000 0.50100000 -137.14070945
-0.00000000 0.50000000 443.79666119
0.00100000 0.49900000 -137.14070945
0.00200000 0.49800000 -359.03904207
0.00300000 0.49700000 359.03903926
0.00400000 0.49600000 137.14071399
0.00500000 0.49500000 -443.79666119
0.00600000 0.49400000 137.14070490
0.00700000 0.49300000 359.03904488
0.00800000 0.49200000 -359.03903645
0.00900000 0.49100000 -137.14071854
0.01000000 0.49000000 443.79666119
0.01100000 0.48900000 -137.14070035
0.01200000 0.48800000 -359.03904769
0.01300000 0.48700000 359.03903364
0.01400000 0.48600000 137.14072308
0.01500000 0.48500000 -443.79666119

補足情報（FW/ツールのバージョンなど）

Visual Studio Code (Version 1.46)を使用しています。

Comments

[maai]

フーリエ級数展開の定義とソースコード中の変数の意味を書いてください（特にMM）
そもそもフーリエ級数展開の求め方が間違っているのでは…？

[ty59rhgaeraoih]

追記いたしました。①~⑤の文が穴埋めになっている形式なのでそれ以外の部分は間違っていません。宜しくお願い致します。

[Penpen7]

本質的な部分ではありませんが、math.hをインクルードしているため、#define M_PI 3.14159265 の定義は必要ありません。
各数学定数を使用するために、
#define _USE_MATH_DEFINES
#include <math.h>
とする必要はあるかもしれませんが。
double fcoeff_bk, fcoeff_akのiはループ変数なのでdoubleではなくintでは？

[Penpen7]

f(t) = 1/2(1 - 2t/T) -2/T < t < T/2
　　 1/2 t = ± T/2
と
f(t) = 1/2{ 1 + cos(2πt / T) -T/2 ≦ t ≦ T/2
の2種類の関数をフーリエ級数展開したいという意味ですかね？(3つ並んでいてわかりづらかったです)
else return 0.5 - t/Tp; //ここは間違っている可能性あり①
まではあってると思いますが、
return 0.5 + 0.5 * cos((2.0 * M_PI * t) / Tp); //ここは間違っている可能性あり②
は別の関数だからいらないのでは？(そもそもこのreturn文は実行されることはなく不要です)

[Penpen7]

学校の課題でしょうから、私からはヒントだけお出しします。
ak = ak + sin(t) / t ;とbk = bk + sin(t) / t ;は間違っています。
ここでは, フーリエ級数展開の係数を求めますが、係数が積分によって定義されているために、数値積分を行います。
数値積分には、区分求積、シンプソンの式、台形公式等がありますが、手っ取り早く実装できるのは区分求積です。ご質問のコードで微小区間の幅をかけていないということは、積分を単なる足し算だと考えていませんか？
もう一度、「区分求積」とフーリエ級数展開の「係数」の定義をみてください。
https://gihyo.jp/dev/serial/01/java-calculation/0072
https://mathtrain.jp/fourierseries

[ty59rhgaeraoih]

ご指導ありがとうございます。
おっしゃる通り、三角波と正弦波の2種類の関数をフーリエ級数展開するプログラムです。
t = -Tp/2.0 + (double)i * dt ;　が何をしているのか分かりません…

[Penpen7]

t = -Tp/2.0 + (double)i * dt ;は
i=0の時, t=-Tp/2.0
i=NNの時, t=-Tp/2.0+NN*dt=Tp/2.0 (∵dt=Tp/NN)
となります。すなわち、-Tp/2.0<=t<Tp/2.0の範囲を最初-Tp/2.0からdtずつ増やしていきながら、いずれTp/2.0にたどり着くと考えてください。(Tp/2.0の方に=はつかないので完全にTp/2.0にはなりませんが簡単のため)
なぜこのようなことをしているのかというと、区分求積では関数f(t)をt=aからt=bまで積分するとすると、f(t)*dt=(関数の値)*(幅)をtをaからbまでdt(微小幅)ずつ「ずらしながら」足し合わせていく操作をするからです。
関数の値を求めるためにはtの値が必要のため、tの値をその文で計算しています。

[ty59rhgaeraoih]

回答ありがとうございます。遅れてしまい申し訳ございません。あれからずっと考えていたのですが、
矩形近似のx_i-x_i-1(横)に当たる部分はdtなのでしょうか。また、y_i-1(縦)にあたる部分はなんでしょうか。
そう考えて　ak = ak+ak_dt(i,t) * dt ;　としてみましたが、うまくいきませんでした。
また、main関数内の
for (k=1; k<=MM; k++) {
aa[k] = fcoeff_ak(k) ; bb[k] = fcoeff_bk(k) ; /* フーリエ展開係数 a_k, b_k を求める */ ;
}
の部分でなぜfor文を使っているのか分かりません。
回答宜しくお願い致します。

[Penpen7]

x_i-x_i-1が積分の方向(dx)を言っているのだとすれば長方形の横幅にあたります。このソースコードであればdtにあたりますね。
y_{i-1}は被積分関数に対応します。先日コメントした参考URL(https://mathtrain.jp/fourierseries)のa_nの式であれば、被積分関数は2/T*f(x)*cos(2*pi*n*x/T)ですのでyにはこれを当てはめることになります。
for (k=1; k<=MM; k++) {
aa[k] = fcoeff_ak(k) ; bb[k] = fcoeff_bk(k) ;
}
は全てのkに対して、結果を得たいためにfor文を使っています。
(Tp/kが周期に対応)

[Penpen7]

ak = ak+ak_dt(i,t) * dt ;であっていると思います。
同様にbk = bk + bk_dt(i,t) * dt;としましたか？

[ty59rhgaeraoih]

回答ありがとうございます。bkの部分も変更しました。コンパイルでエラーは出ませんが、実行しても何も出力されません…
それと、fxfreqと名前を付けた関数内の
for ( k = 0 ; k <= MM ; k++ ) ff = *aa*cos(((double)k * 2.0 * M_PI * t)/ Tp ) + *bb*sin(((double)k * 2.0 * M_PI * t)/Tp ) ;
は正しいですか？
回答宜しくお願い致します。

[Penpen7]

ソースコードを質問の編集でもう一度貼り直していただけませんか？
自分の手で動かして確認したいので...
fxfreqは何をする部分でしょうか？おそらく求めたフーリエ係数からf(t)を再現する作業をしているのですかね？
(fhとftが同じになるはずとのことからそのように予想できますが)

[Penpen7]

もしその意味なら*aaはaa[k], *bbはbb[k]とするのが正しいです。
配列を引数に取りたいときはポインタを使うことがあります。

[ty59rhgaeraoih]

かしこまりました。よろしくお願いします。
fxfreqはフーリエ展開係数 a_k,b_kを求めています。

[ty59rhgaeraoih]

#include <stdio.h>
#include <math.h>
#define M_PI 3.14159265
#define MM 300
#define Tp 1.0

double aa[MM+1], bb[MM+1] ;
int NN = 3 * MM ;

double ft(double t)
{
double zero = 1.0e-8 ;
while ( t < -Tp/2.0 ){ t = t+Tp ;if( t > -Tp/2.0 ) break ;}
while ( t > Tp/2.0 ){ t = t-Tp ;if( t < Tp/2.0 ) break ;}

if( fabs( t+Tp/2.0 ) < zero || fabs( t-Tp/2.0) < zero ) return 0.5 ;
else return 0.5 - t/Tp;

return 0.5 + 0.5 * cos((2.0 * M_PI * t) / Tp);
}

double ak_dt(int k, double t)
{
return ft(t) * cos( (double)k * 2.0 * M_PI / Tp * t ) ;
}

double fcoeff_ak( int k )
{
double t, i, ak = 0.0, dt = Tp/(double)NN ;
for( i = 0 ; i < NN ; i ++ ) {
t = -Tp/2.0 + (double)i * dt ;
ak = ak + ak_dt(i,t) * dt ;
}

return ak ;
}

double bk_dt(int k, double t)
{ /* フーリエ展開係数 b_k を求めるときの被関数関数 */
return ft(t) * sin( (double)k * 2.0 * M_PI / Tp * t ) ;
}

double fcoeff_bk(int k)
{ /* 関数 f(t) のフーリエ展開係数 b_k を数値積分で求める */
double t, i, bk = 0.0, dt = Tp/(double)NN ; // 初期化 bk=0.0 を忘れずに!
for( i = 0 ; i < NN ; i ++ ) {
t = -Tp/2.0 + (double)i * t ;
bk = bk + bk_dt(i,t) * dt ; // 矩形近似，または台形近似による積分
}
return bk ;
}

double fxfreq( double *aa, double *bb, int M, double t)
{ /* M 次までのフーリエ展開係数を用いて，ある t 点での f(t) を計算する */
int k ; double ff = 0.0 ;
for ( k = 0 ; k <= MM ; k++ ) ff = aa[k]*cos(((double)k * 2.0 * M_PI * t)/ Tp ) + bb[k]*sin(((double)k * 2.0 * M_PI * t)/Tp ) ;
return ff;
}

int main()
{
int i, k, n=2 ; double t, fh ;
int KK = 1000 ; /* 1周期分のグラフ描写の点数 */

aa[0] = fcoeff_ak(0)/2.0 ; bb[0] = 0.0 ; /* フーリエ展開係数 a_0, b_0 を求める */ ;
for (k=1; k<=MM; k++) {
aa[k] = fcoeff_ak(k) ; bb[k] = fcoeff_bk(k) ; /* フーリエ展開係数 a_k, b_k を求める */ ;
}

/* 求めたフーリエ展開係数から n 周期分の波形 fh を再構成・描写する */
for (i=0; i<= n*KK; i++) {
t = -Tp*(double)n/2.0 + Tp/(double)KK * (double)i ;
fh = fxfreq( aa, bb, MM, t ) ;
printf(" %10.8f %10.8f %10.8f\n", t, ft(t), fh );
}
return 0;
}

Answer 1

引数を間違えている箇所等がありました。(見比べてください)
デバッグのため、ソースコードを少し改変しています。

#include <math.h>
#include <stdio.h>
#define MM 100

const double Tp = 1.0;
const int NN = 3*MM;
const double dt = Tp / (double)NN;
double aa[MM + 1], bb[MM + 1];

double ft(double t)
{
  double zero = 1.0e-8;
  while (t < -Tp / 2.0) {
    t = t + Tp;
    if (t > -Tp / 2.0)
      break;
  }
  while (t > Tp / 2.0) {
    t = t - Tp;
    if (t < Tp / 2.0)
      break;
  }

  if (fabs(t + Tp / 2.0) < zero || fabs(t - Tp / 2.0) < zero)
    return 0.5;
  else
    return 0.5 - t / Tp;
}

double ak_dt(int k, double t)
{
  return ft(t) * cos((double)k * 2.0 * M_PI / Tp * t);
}

double fcoeff_ak(int k)
{
  double ak = 0.0;
  for (int i = 0; i < NN; i++) {
    double t = -Tp / 2.0 + (double)i * dt;
    ak += ak_dt(k, t) * dt;
  }
  ak *= 2 / Tp;
  return ak;
}

double bk_dt(int k, double t)
{ /* フーリエ展開係数 b_k を求めるときの被関数関数 */
  return ft(t) * sin((double)k * 2.0 * M_PI / Tp * t);
}

double fcoeff_bk(int k)
{ /* 関数 f(t) のフーリエ展開係数 b_k を数値積分で求める */
  double bk = 0.0;
  for (int i = 0; i < NN; i++) {
    double t = -Tp / 2.0 + (double)i * dt;
    bk += bk_dt(k, t) * dt; // 矩形近似，または台形近似による積分
  }
  bk *= 2 / Tp;
  return bk;
}

double fxfreq(double* aa, double* bb, int M, double t)
{ /* M 次までのフーリエ展開係数を用いて，ある t 点での f(t) を計算する */
  double ff = 0.0;
  for (int k = 0; k <= MM; k++)
    ff += aa[k] * cos(((double)k * 2.0 * M_PI * t) / Tp) + bb[k] * sin(((double)k * 2.0 * M_PI * t) / Tp);
  return ff;
}

int main()
{
  const int KK = 1000; /* 1周期分のグラフ描写の点数 */

  aa[0] = fcoeff_ak(0) / 2.0;
  bb[0] = 0.0; /* フーリエ展開係数 a_0, b_0 を求める */
  for (int k = 1; k <= MM; k++) {
    aa[k] = fcoeff_ak(k);
    bb[k] = fcoeff_bk(k); /* フーリエ展開係数 a_k, b_k を求める */
  }

  /* 求めたフーリエ展開係数から n 周期分の波形 fh を再構成・描写する */
  for (int i = 0; i <= KK; i++) {
    double t = -Tp * 0.5 + Tp / (double)KK * (double)i;
    double fh = fxfreq(aa, bb, MM, t);
    printf(" %10.8f %10.8f %10.8f\n", t, ft(t), fh);
  }
  return 0;
}