8回目　課題　解説

#include <stdio.h>
#include <math.h>

struct Complex{
　double re; /* 実部 */
　double im; /* 虚部 */
};


void clear_array(double x[], int n); // 要素数 n の配列 x に，すべて 0 を代入する（初期化）．
void square_wave(double x[], int n); // 要素数 n の配列 x に，任意の周波数（例えば 5 ）の矩形波を書き込む．振幅は任意，例えば -1 から 1 など．
void sin_waves(double x[], int n, double f1, double f2, double f3); // 配列に周波数 f1, f2, f3，計3個の正弦波または余弦波の和を入れる．
Complex ToComplex(double x, double y);
Complex Cadd(Complex a, Complex b);
Complex Cmul(Complex a, Complex b);
Complex Cexp(Complex a);
double Cabs(Complex a);
void dft(const double in[], Complex out[], const int N);


int main(){
　const int n = 100;
　double x1[n] = {0}, x2[n] = {0}, x11[n] = {0}, x21[n] = {0};
　Complex dx1[n], dx2[n];
　FILE *fp;

　//配列の初期化
　clear_array(x1, n);
　clear_array(x2, n);
　clear_array(x11,n);
　clear_array(x21,n);

　//配列に各波形の値を代入
　square_wave(x1, n);
　sin_waves(x2, n, 1, 2, 5);

　//dftの計算
　dft(x1, dx1, n);
　dft(x2, dx2, n);

　//スペクトルの計算
　for(int i = 0; i < n; i++){
　　x11[i] = Cabs(dx1[i]) * Cabs(dx1[i]);
　　x21[i] = Cabs(dx2[i]) * Cabs(dx2[i]);
　}

　//ファイル出力
　fp = fopen("153R000000-08-1.csv", "w");
　if(fp == NULL){
　　printf("file cannot be opened\n");
　　return -1;
　}

　fprintf(fp, "周波数,矩形波,合成波\n");
　for(int i = 0; i < n/2; i++){
　　fprintf(fp, "%d,%f,%f\n", i, x11[i], x21[i]);
　}

　fclose(fp);
　printf("ファイル出力完了!\n");

　return 0;
}

void clear_array(double x[], int n){
　for(int i = 0; i < n; i++){
　　x[i] = 0;
　}
}

void square_wave(double x[], int n){
　int freq = 10; //任意の半周波数をここで設定
　for(int i = 0; i < n; i++){
　　if((i/freq)%2 == 0) x[i] = -1; //iを任意の半周期で割った解の整数の位が2の倍数なら-1を代入,それ以外で1を代入
　　else x[i] = 1;
　}
}

void sin_waves(double x[], int n, double f1, double f2, double f3){
　for(int i = 0; i < n; i++){
　　x[i] = (2.0 * sin(2*M_PI*i/n * f1) //M_PIはmath.hに定義される円周率の定数
　　　+ 0.5 * cos(2*M_PI*i/n * f2)
　　　+ 0.25 * sin(2*M_PI*i/n * f3));
　}
}

Complex ToComplex(double x, double y){ /*複素数のﾊﾟﾗﾒｰﾀを与える*/
　Complex z; /*Complex型の変数を用意*/
　z.re = x;
　z.im = y;
　return z; /*構造体として返す（実部，虚部をパラメータとして持つ）*/
}

Complex Cadd(Complex a, Complex b){ /*複素数の足し算*/
　Complex z;
　z.re = a.re + b.re;/*複素数aの実部と複素数bの実部を足す．*/
　z.im = a.im + b.im;
　return z;
}

Complex Cmul(Complex a, Complex b){ /*複素数の掛け算*/
　Complex z;
　z.re = (a.re * b.re) - (a.im * b.im);
　z.im = (a.re * b.im) + (b.re * a.im);
　return z;
}

Complex Cexp(Complex a){ /*指数関数*/
　Complex z;
　z.re = exp(a.re) * cos(a.im);
　z.im = exp(a.re) * sin(a.im);
　return z;
}

double Cabs(Complex a){ /*絶対値*/
　double x;
　x = sqrt(a.re * a.re + a.im * a.im);
　return x;
}

void dft(const double in[], Complex out[], const int N){
//
// 離散フーリエ変換(DFT)
// 元データは実数，変換後は複素数
//
// in : 元データ．これは変更しない．実数の配列．
// out : 計算結果．こちらは複素数．
// N : データ個数．inもoutも同じデータ数でなければならない
//
　int t;
　/* 出力用の配列を初期化 */
　for(t=0; t<N; t++) {
　　out[t].re = 0.0;
　　out[t].im = 0.0;
　}
　for(int f=0; f<N; f++) {
　// 以下のループで 1個のF[omega]を計算する．
　　double omega = 2 * M_PI * f;
　　for(t=0; t<N; t++) {
　　　/* 積分の中身の計算 */
　　　out[f] = Cadd(out[f], Cmul( ToComplex(in[t], 0.0), Cexp( ToComplex(0.0, -omega * t/N) )));
　　}
　}
}


構造体





関数の宣言










メイン関数












































配列初期化





矩形波の関数






正弦波合成の関数


#include <stdio.h>
#include <math.h>

struct Complex{
　double re; /* 実部 */
　double im; /* 虚部 */
};


void clear_array(double x[], int n); // 要素数 n の配列 x に，すべて 0 を代入する（初期化）．
Complex ToComplex(double x, double y);
Complex Cadd(Complex a, Complex b);
Complex Cmul(Complex a, Complex b);
Complex Cexp(Complex a);
double Cabs(Complex a);
void dft(const double in[], Complex out[], const int N);



int main(){
　const int N = 1000; //要素数
　double data[N] = {0}, spdata[N] = {0}; //dataは元データ、spdataはパワースペクトルのデータ
　Complex dftdata[N] = {0}; //dftでの計算結果を格納するComplex型配列
　FILE *ifp, *ofp; //FILE構造体へのポインタ,ifp:input, ofp:output
　char ifilename[] = "wave.txt", ofilename[] = "153R000000-08-2.csv"; /* ファイル名 */

　ifp = fopen(ifilename, "r"); /* 読み込みモードで開く */
　if(ifp == NULL){
　　printf("入力ファイルが開けない\n");
　　return -1; /* エラー終了 */
　}

　for(int i=0; i<N; i++){
　　fscanf(ifp, "%lf", &data[i]);
　}
　fclose(ifp);

　/* 念のため画面に出力して確認 */
　for(int i=0; i<N; i++){
　　printf("data[%d] = %lf\n", i, data[i]);
　}

　//spdataを初期化
　clear_array(spdata, N);

　//dft演算
　dft(data, dftdata, N);

　//パワースペクトルを計算
　for(int i = 0; i < N; i++){
　　spdata[i] = Cabs(dftdata[i]) * Cabs(dftdata[i]);
　}

　//ファイル出力
　ofp = fopen(ofilename, "w");
　fprintf(ofp, "周波数,パワースペクトル\n");
　if(ofp == NULL){
　　printf("出力ファイルが開けない\n");
　　return -1;
　}

　for(int i = 0; i < N/2; i++){
　　fprintf(ofp, "%d,%f\n", i, spdata[i]);
　}

　fclose(ofp);
　printf("ファイル出力完了\n");

　return 0;
}

void clear_array(double x[], int n){
　for(int i = 0; i < n; i++){
　　x[i] = 0;
　}
}

Complex ToComplex(double x, double y){ /*複素数のﾊﾟﾗﾒｰﾀを与える*/
　Complex z; /*Complex型の変数を用意*/
　z.re = x;
　z.im = y;
　return z; /*構造体として返す（実部，虚部をパラメータとして持つ）*/
}

Complex Cadd(Complex a, Complex b){ /*複素数の足し算*/
　Complex z;
　z.re = a.re + b.re;/*複素数aの実部と複素数bの実部を足す．*/
　z.im = a.im + b.im;
　return z;
}

Complex Cmul(Complex a, Complex b){ /*複素数の掛け算*/
　Complex z;
　z.re = (a.re * b.re) - (a.im * b.im);
　z.im = (a.re * b.im) + (b.re * a.im);
　return z;
}

Complex Cexp(Complex a){ /*指数関数*/
　Complex z;
　z.re = pow(exp(1), a.re) * cos(a.im);
　z.im = pow(exp(1), a.re) * sin(a.im);
　return z;
}

double Cabs(Complex a){ /*絶対値*/
　double x;
　x = sqrt(a.re * a.re + a.im * a.im);
　return x;
}

void dft(const double in[], Complex out[], const int N){
//
// 離散フーリエ変換(DFT)
// 元データは実数，変換後は複素数
//
// in : 元データ．これは変更しない．実数の配列．
// out : 計算結果．こちらは複素数．
// N : データ個数．inもoutも同じデータ数でなければならない
//
　int t;
　/* 出力用の配列を初期化 */
　for(t=0; t<N; t++) {
　　out[t].re = 0.0;
　　out[t].im = 0.0;
　}
　for(int f=0; f<N; f++) {
　// 以下のループで 1個のF[omega]を計算する．
　　double omega = 2 * M_PI * f;
　　for(t=0; t<N; t++) {
　　　/* 積分の中身の計算 */
　　　out[f] = Cadd(out[f], Cmul( ToComplex(in[t], 0.0), Cexp( ToComplex(0.0, -omega * t/N) )));
　　}
　}
}


構造体





関数の宣言









メイン関数



















































配列初期化