Спектр периодической последовательности прямоугольных импульсов

#include <stdio.h>
#include <string.h>
#include "dspl.h"

/* Размер векторов входных сигналов и огибающей спетра */
#define N  1000 

/* Период повторения импульса. Для изменения скважности мы будем менять 
 * длительность импульса при фиксированном периоде повторения */
#define T  4.0

/* Амплитуда */
#define A  2.0

/* Количество спектральных гармоник разложения в ряд Фурье  */
#define M  41


/* длина команды Gnuplot */
#define PLOTCMD_LEN 256


int main(int argc, char* argv[])
{

  double    t1[N];   /* время (сек) на одном периоде повторения      */
  double    t4[N];   /* время (сек) на четырех периодах повторения   */
  double    s[N];    /* входной сигнал                               */
  complex_t S[M];    /* комплексный спектр периодического сигнала    */
  double    Smag[M]; /* амплитудный спектр периодического сигнала    */
  double    w[M];    /* частота (рад/c) дискретного спектра          */
  double    wc[N];   /* частота (рад/с) огибающей спектра            */
  double    Sc[N];   /* огибающая спектра                            */
  double    tau;     /* длительность импульса                        */

  /* скважность */
  double    Q[3] = {5.0, 2.0, 1.25};

  int q, m, n;
  
  char fname[64];            /* имя файла данных */
  char plotcmd[PLOTCMD_LEN]; /* Команда Gnuplot  */ 
  
  void* hdspl;  /* DSPL handle        */
  void* hplot;  /* GNUPLOT handle     */
  
  hdspl = dspl_load();
  if(!hdspl)
  {
    printf("Cannot to load libdspl!\n");
    return 0;
  }
  
  /* Вектор частот непрерывной огибаюхей вида sin(w/2*tau) / (w/2*T) */
  linspace(-M_PI*(double)M/(double)T, M_PI*(double)M/(double)T,  N, DSPL_SYMMETRIC, wc);
  
  /* заполнение массива временных отсчетов  */
  /* на одном периоде повторения сигнала    */
  linspace(-T/2.0, T/2.0,  N, DSPL_PERIODIC, t1);
  
  /* заполнение массива временных отсчетов
   * на 4-x периодах повторения сигнала
   * для отображения на осциллограмме */
  linspace(-T*2.0, T*2.0,  N, DSPL_PERIODIC, t4);
  
  
  /* Построение графиков пакетом GNUPLOT */
  gnuplot_create(argc, argv, 800, 640,  "img/fourier_series_rec.png", &hplot);
  gnuplot_cmd(hplot, "unset key");
  gnuplot_cmd(hplot, "set multiplot layout 3,2 rowsfirst");
  gnuplot_cmd(hplot, "set yrange [0:2.2]");
  
  for(q = 0; q < 3; q++)
  {
    tau = T/Q[q];

    /* 4 периода повторения п-импульса скважности Q[q] */
    signal_pimp(t4, N, A, tau, 0.0, T, s);

    /* сохранение в текстовый файл временных осциллограмм */
    sprintf(fname, "dat/pimp_time_%.2lf.csv", Q[q]);
    writetxt(t4, s, N, fname);

    /* Построение временнОй осциллограммы */
    sprintf(plotcmd, "plot '%s' with lines", fname);
    gnuplot_cmd(hplot, plotcmd);

    /* один период повторения п-импульса скважности Q[q] */
    signal_pimp(t1, N, A, tau, 0.0, T, s);

    /* разложение в ряд Фурье */
    fourier_series_dec(t1, s, N, T, M, w, S);

    /* Рассчет амплитудного спектра */
    for(m = 0; m < M; m++)
    {
      /*printf("S[%d] = %f     %f\n", m, RE(S[m]), IM(S[m]));*/
      Smag[m] = ABS(S[m]);
    }

    /* Сохранение в файл амплитудного спетра для скважности Q[q] */
    sprintf(fname, "dat/pimp_freq_discrete_%.2lf.csv", Q[q]);
    writetxt(w, Smag, M, fname);

    /* Построение на график амплитудного спектра для заданной скважности */
    sprintf(plotcmd, "plot '%s' with  impulses lt 1 ,\\", fname);
    printf("%s\n", plotcmd);
    gnuplot_cmd(hplot, plotcmd);

    sprintf(plotcmd, "'%s'  with points pt 7 ps 0.5 lt 1 ,\\", fname);
    printf("%s\n", plotcmd);
    gnuplot_cmd(hplot, plotcmd);

    /* Расчет огибающей */
    for(n = 0; n < N; n++)
      Sc[n] = (wc[n] == 0.0) ? A/Q[q] : fabs( A * sin(0.5*wc[n]*tau) / (0.5*wc[n] * T));

    /* сохранение огибающей в файл для скважности Q[q] */
    sprintf(fname, "dat/pimp_freq_cont_%.2lf.csv", Q[q]);
    writetxt(wc, Sc, N, fname);
    
    /* Построение на график непрерывной огибающей 
    амплитудного спектра для заданной скважности */
    sprintf(plotcmd, "'%s'  with lines", fname);
    printf("%s\n", plotcmd);
    gnuplot_cmd(hplot, plotcmd);
  }

  gnuplot_cmd(hplot, "unset multiplot");
  gnuplot_close(hplot);
  
  /* remember to free the resource */
  dspl_free(hdspl);
  return 0;
}