/*Maximaでの例*/

/*maximaの命令（関数）や記述方法を覚える必要はないが、目的をこの程度で表せるという点には注目*/

/*標本化定理、本当は、例：http://www.image.med.osaka-u.ac.jp/member/yoshi/ouec_lecture/digital_processing/handout/Sampling_theorem.pdf*/

/*●標本化定理の理解深化のため*/
/*電話と同じ8000Hzでサンプリングしたとして、それを密な棒グラフ的に再生したときの波形を解析してみる*/
/*するろ分かるが、それは、元のアナログ信号と同じ周波数以外に、サンプリング周波数の半分以上の周波数に余計な高周波の波形が含まれている*/

/*一旦、maximaに定義したものを全てリセット*/
kill(all);

/*プロット用の関数scatterplotを使うためにパッケージを読み込む*/
load(descriptive); 

/*元の信号の例：適当な振幅の、500Hzと1000Hzと2000Hzの信号からなっている場合*/
s(t):=sin(2*%pi*500*t)+1/2*sin(2*%pi*1000*t)+1/4*sin(2*%pi*2000*t);

/*それぞれの波形と合わせた波形を描画*/
plot2d( [sin(2*%pi*500*t),1/2*sin(2*%pi*1000*t),1/4*sin(2*%pi*2000*t),s(t)], [t,0,1/500]);
/*合わせた波形だけを描画*/
plot2d( [s(t)], [t,0,1/500]);

/*sin(2*%pi*fx*t)の成分がどれだけ量あるか（フーリエ係数）を調べる。かけて1秒間分積分すればよい*/
/*代わりに、適当な間隔（ここでは1/8000秒）ごとの面積の足し算を1/50秒間分計算したのを使って、求めている。少し時間がかかる*/
/*「2*」は数学的に2倍して係数になるから。「*50」は1/50秒分のを1秒分にするため。*/
/*floatは数式を実際に計算させる関数*/
/*グラフは横軸が周波数、縦軸がフーリエ係数。グラフのピークは実は幅がなく、非ゼロのプロットで描画している*/
scatterplot(
  makelist( [fx,float(2*integrate( s(t)*sin(2*%pi*fx*t),t,0,1/50)*50)] , fx,250,8000,50)
  ,point_type = diamant
);　

/*描画範囲を広げてみた*/
scatterplot(
  makelist( [fx,float(2*integrate( s(t)*sin(2*%pi*fx*t),t,0,1/50)*50)] , fx,250,24000,250)
  ,point_type = diamant
);

/*8000Hzで標本化して（ここでは標本化だけに注目するので、量子化誤差が無視できるとして）符号化して、その情報を使って再生した場合の信号の式*/
/*8000分の1秒ごとに値が変わるような関数にすればよい。そのために、8000分の1秒間隔で切り捨てにした時間で関数を使っている*/
d(t):=s( ( floor( t / (1/8000) ) * (1/8000) ) );
/**/ 
plot2d( [s(t),d(t)], [t,0,1/500]);

/*以下の数行は、1/8000秒間での平均値を符号化して復号でdとして用いることができた場合の解析用の関数（再）定義*/
xd(t):=integrate( s(x), x, ( floor( t / (1/8000) ) * (1/8000) ), ( floor( t / (1/8000) ) * (1/8000) ) + (1/8000) ) *8000;
d2all(t):=( ( 1/(2*%pi*500)*( -cos(2*%pi*500*(t+1/8000)) +cos(2*%pi*500*(t)) ) )*1 +
            ( 1/(2*%pi*1000)*( -cos(2*%pi*1000*(t+1/8000)) +cos(2*%pi*1000*(t)) ) )*(1/2) +
            ( 1/(2*%pi*2000)*( -cos(2*%pi*2000*(t+1/8000)) +cos(2*%pi*2000*(t)) ) )*(1/4) 
          )    *8000;
d(t):=d2all( ( floor( t / (1/8000) ) * (1/8000) ) );
plot2d( [s(t),d(t)], [t,0,1/500]);


/*sin(2*%pi*fx*t)の成分がどれだけあるかを調べる。かけて1秒間分積分すればよい*/
/*さっきの式で、sをdに変えるだけでは、maximaでは計算されなかったので、一工夫*/
/*sをdに変えて、1/8000秒ごとの積分を1/50秒分加算する。*/
/*スピードアップのために、dの項は積分の外に出して、sinは積分関数-cos(2*%pi*fx*t)を使って求めた*/
scatterplot(
  makelist( [fx,
             float(2*sum( float( 
               d(i*1/8000)* 
               1/(2*%pi*fx)*( -cos(2*%pi*fx*((i+1)*1/8000)) +cos(2*%pi*fx*(i*1/8000)) )   
             ),i,0,8000/50-1)*50)
            ], fx,250, 24000,250)
  ,point_type = diamant
);

/* 変形途中の式 makelist( [fx,float(2*sum( float( d(i*1/8000)* integrate(sin(2*%pi*fx*t),t,i*1/8000,(i+1)*1/8000) ) ,i,0,8000/50-1)*50/8000)] , fx,1000, 1000,250) */

/*正確には、sin(2*%pi*fx*t)の成分とcos(2*%pi*fx*t)の成分があることを考慮する。少し時間がかかる*/
/*若干値が異なり、高周波の振幅が小さくなるとすれば、厳密ばくし歯の関数を使っていないから*/
scatterplot(
  makelist( [fx,
             (
              (
               float(2*sum( float( 
                 d(i*1/8000)* 
                 1/(2*%pi*fx)*( -cos(2*%pi*fx*((i+1)*1/8000)) +cos(2*%pi*fx*(i*1/8000)) )   
               ),i,0,8000/50-1)*50)
              )^2+
              (
               float(2*sum( float( 
                 d(i*1/8000)* 
                 1/(2*%pi*fx)*(  sin(2*%pi*fx*((i+1)*1/8000)) -sin(2*%pi*fx*(i*1/8000)) )   
               ),i,0,8000/50-1)*50)
              )^2
             )^0.5
            ], fx,250, 24000,250)
  ,point_type = diamant
);


/*以下の数行は、計算はしているが、デルタ関数のようになるためか、0以外の値がプロットされなかった*/
/*fに対するd(t)の振幅y2(f)を、cos成分、sin成分より求める。但し、参考用*/
/*y(f):=( coeff(d(t),cos(2*%pi*f*t),1)^2 + coeff(d(t),sin(2*%pi*f*t),1)^2 )^0.5;*/
/*fの適当な範囲についてy(f)を描画。但し、上記の定数では約1MHzで同調することが既知であるのを利用*/
/*y(f):=2*( coeff(sin(2*%pi*440*t),cos(2*%pi*f*t),1)^2 + coeff(sin(2*%pi*440*t),sin(2*%pi*f*t),1)^2 )^0.5*/
/*plot2d( [y(f)], [f,439.9*10^0,440.1*10^0])*/

/**/

ーーーーー
/*●標本化定理の理解深化のため*/
/*電話と同じ8000Hzでサンプリングしたとして、それを密な棒グラフ的に再生したときの波形を解析してみる*/
/*それは、元のアナログ信号と同じ周波数以外に、サンプリング周波数の半分以上の周波数に余計な高周波の波形が含まれている*/

/*一旦、maximaに定義したものを全てリセット*/
kill(all);

/*プロット用の関数scatterplotを使うためにパッケージを読み込む*/
load(descriptive); 

/*元の信号の例：適当な振幅の、500Hzと1000Hzと2000Hzの信号からなっている場合*/
s(t):=sin(2*%pi*500*t)+1/2*sin(2*%pi*1000*t)+1/4*sin(2*%pi*2000*t);

/*それぞれの波形と合わせた波形を描画*/
plot2d( [sin(2*%pi*500*t),1/2*sin(2*%pi*1000*t),1/4*sin(2*%pi*2000*t),s(t)], [t,0,1/500]);
/*合わせた波形だけを描画*/
plot2d( [s(t)], [t,0,1/500]);

/*sin(2*%pi*fx*t)の成分がどれだけ量あるか（フーリエ係数）を調べる。かけて1秒間分積分すればよい*/
/*代わりに、適当な間隔（ここでは1/8000秒）ごとの面積の足し算を1/50秒間分計算したのを使って、求めている。少し時間がかかる*/
/*「2*」は数学的に2倍して係数になるから、「*50」は1/50秒分のを1秒分にするため*/
/*グラフは横軸が周波数、縦軸がフーリエ係数。グラフのピークは実は幅がなく、非ゼロのプロットで描画している*/
scatterplot(
  makelist( [fx,float(2*integrate( s(t)*sin(2*%pi*fx*t),t,0,1/50)*50)] , fx,250,8000,50)
  ,point_type = diamant
);　

/*描画範囲を広げてみた*/
scatterplot(
  makelist( [fx,float(2*integrate( s(t)*sin(2*%pi*fx*t),t,0,1/50)*50)] , fx,250,24000,250)
  ,point_type = diamant
);

/*8000Hzで標本化して（ここでは標本化だけに注目するので、量子化誤差が無視できるとして）符号化して、その情報を使って再生した場合の信号の式*/
/*8000分の1秒ごとに値が変わるような関数にすればよい。そのために、8000分の1秒間隔で切り捨てにした時間で関数を使っている*/
d(t):=s( ( floor( t / (1/8000) ) * (1/8000) ) );
/**/ 
plot2d( [s(t),d(t)], [t,0,1/500]);

/*以下の数行は、1/8000秒間での平均値を符号化して復号でdとして用いることができた場合の解析用の関数（再）定義*/
xxxd(t):=integrate( s(x), x, ( floor( t / (1/8000) ) * (1/8000) ), ( floor( t / (1/8000) ) * (1/8000) ) + (1/8000) ) *8000;
d2all(t):=( ( 1/(2*%pi*500)*( -cos(2*%pi*500*(t+1/8000)) +cos(2*%pi*500*(t)) ) )*1 +
            ( 1/(2*%pi*1000)*( -cos(2*%pi*1000*(t+1/8000)) +cos(2*%pi*1000*(t)) ) )*(1/2) +
            ( 1/(2*%pi*2000)*( -cos(2*%pi*2000*(t+1/8000)) +cos(2*%pi*2000*(t)) ) )*(1/4) 
          )    *8000;
ad(t):=d2all( ( floor( t / (1/8000) ) * (1/8000) ) );
plot2d( [s(t),d(t)], [t,0,1/500]);


/*sin(2*%pi*fx*t)の成分がどれだけあるかを調べる。かけて1秒間分積分すればよい*/
/*さっきの式で、sをdに変えるだけでは、maximaでは計算されなかったので、一工夫*/
/*sをdに変えて、1/8000秒ごとの積分を1/50秒分加算する。*/
/*スピードアップのために、dの項は積分の外に出して、sinは積分関数-cos(2*%pi*fx*t)を使って求めた*/

/*正確には、sin(2*%pi*fx*t)の成分とcos(2*%pi*fx*t)の成分があることを考慮する。少し時間がかかる*/
/*若干値が異なり、高周波の振幅が小さくなるとすれば、厳密ばくし歯の関数を使っていないから*/


/*dを1/8000の更に1/mの幅だけm倍の振幅の関数として用いたとき*/
m:4;
scatterplot(
  makelist( [fx,
             (
              (
               float(2*sum( float( 
                 m*d(i*1/8000)* 
                 1/(2*%pi*fx)*( -cos(2*%pi*fx*((i+1/m)*1/8000)) +cos(2*%pi*fx*(i*1/8000)) )   
               ),i,0,8000/50-1)*50)
              )^2+ 
              (
               float(2*sum( float( 
                 m*d(i*1/8000)* 
                 1/(2*%pi*fx)*(  sin(2*%pi*fx*((i+1/m)*1/8000)) -sin(2*%pi*fx*(i*1/8000)) )   
               ),i,0,8000/50-1)*50)
              )^2
             )^0.5
            ], fx,250, 24000,250)
  ,point_type = diamant
);

/**
平均値を標本値にしたらどう変わるか。
まだやっていない。
/**/


/*●OFDMの理解のため*/


/*一旦、maximaに定義したものを全てリセット*/
kill(all);

/*※どんな計算をしているかに注目。式の入力方法にはとらわれないこと。*/
f0:25*10^8;
b(t):=if t<=1/1000 then 1 else 0;
s0(t):=sin(2*%pi*f0*t)*b(t);

/* error ni
plot2d( integrate(s0(t)*sin(2*%pi*fx*t),t,0,1), [fx,f0-2000,f0+2000]);
plot2d( integrate(s0(t)*sin(2*%pi*fx*t),t,0,1/1000), [fx,f0-2000,f0+2000]);
*/

/* sin 成分のみ */
plot2d( integrate(sin(2*%pi*f0*t)*sin(2*%pi*fx*t),t,0,1/1000), [fx,f0-2000,f0+2000]);

/* cos 成分のみ */
plot2d( integrate(sin(2*%pi*f0*t)*cos(2*%pi*fx*t),t,0,1/1000), [fx,f0-2000,f0+2000]);

/* sin^2 + cos^2 成分 */
plot2d( 
( 
( integrate(sin(2*%pi*f0*t)*sin(2*%pi*fx*t),t,0,1/1000) )^2+
( integrate(sin(2*%pi*f0*t)*cos(2*%pi*fx*t),t,0,1/1000) )^2
)^0.5, [fx,f0-2000,f0+2000]);

/* error ni
plot2d( integrate(s0(t)*sin(2*%pi*fx*t),t,0,1), [fx,f0-2000,f0+2000]);]
*/


/* 飛び飛びの周波数で積分値を求めるような式にした場合 */
load(descriptive); 

/* sin 成分のみ */
scatterplot(
makelist( [fx,float(2*integrate(sin(2*%pi*f0*t)*sin(2*%pi*fx*t),t,0,1/1000))] , fx,f0-2000,f0+2000,250 )
,point_type = diamant
);

/* sin^2 + cos^2 成分 */
scatterplot(
makelist( [fx,float(
(
( 2*integrate(sin(2*%pi*f0*t)*sin(2*%pi*fx*t),t,0,1/1000) )^2+
( 2*integrate(sin(2*%pi*f0*t)*cos(2*%pi*fx*t),t,0,1/1000) )^2
)^0.5
)] , fx,f0-2000,f0+2000,250 )
,point_type = diamant
);

/*OFDM終わり*/


/*●AM変調の理解のため  以下は、数値計算では理論的考察と少し違う結果が出てしまうことの例になっている*/

/*一旦、maximaに定義したものを全てリセット*/
kill(all);

/*※どんな計算をしているかに注目。式の入力方法にはとらわれないこと。*/
/*例として、上のド 523.23Hz 付近の音の波形*/
b(t):=sin(2*%pi*500*t);
/*搬送波と呼ばれる高い周波数の信号とのかけ合わせ 参考：MRTラジオ936kHz*/
f0:1*1000000;
s0(t):=sin(2*%pi*f0*t)*(2+b(t));

plot2d( s0(t), [t,0,1/100]);
plot2d( s0(t), [t,0,1/10000]);

/* sin 成分のみ */
plot2d( integrate( 100* s0(t)*sin(2*%pi*fx*t),t,0,1/1), [fx,f0-1000,f0+1000]);

/* cos 成分のみ */
plot2d( integrate( 100* s0(t)*cos(2*%pi*fx*t),t,0,1/100), [fx,f0-1000,f0+1000]);

/* sin^2 + cos^2 成分 */
plot2d( 
( 
( integrate( 100* s0(t)*sin(2*%pi*fx*t),t,0,1/100) )^2+
( integrate( 100* s0(t)*cos(2*%pi*fx*t),t,0,1/100) )^2
)^0.5, [fx,f0-1000,f0+1000]);


/*以下は積分ができなかった　原因不明　 sin(2*%pi*f0*t)*　を付けると処理されなかった*/

/* 飛び飛びの周波数で積分値を求めるような式にした場合 */
load(descriptive); 

/* sin 成分のみ */
scatterplot(
makelist( [fx,float(100* s0(t) *sin(2*%pi*fx*t),t,0,1/100))] , fx,f0-1000,f0+1000,250 )
,point_type = diamant
);

/* sin^2 + cos^2 成分 */
scatterplot(
makelist( [fx,float(
(
( 2*integrate(s0(t)*sin(2*%pi*fx*t),t,0,1/1000) )^2+
( 2*integrate(s0(t)*cos(2*%pi*fx*t),t,0,1/1000) )^2
)^0.5
)] , fx,f0-2000,f0+2000,250 )
,point_type = diamant
);

/*AM変調終わり*/