locked
FIR and spectrum tools for VC2005 Express RRS feed

  • General discussion

  • afrtohtcoefs -program for converting AFR coeffs to h(t) coeffs for FIR

    // afrtohtcoefs.cpp : Defines the entry point for the console application.
    //

    #include "stdafx.h"
    #include <stdio.h>
    #include <tchar.h>
     #include <iostream>
    #include <stdlib.h>
    #include <cmath>

    #define M_PI 3.1415926535897932384626433832795
    //int _tmain(int argc, _TCHAR* argv[])
    //{
    // return 0;
    //}

     //afrtohtcoefs.cpp


    #define  NMAX 512
    FILE *fp, *fp1; 
     
    int k,n;
    int w,t,N;
    double h[NMAX];
    double h_im[NMAX];
    double Re_H[NMAX];
    double Im_H[NMAX];
    double Re_x[NMAX];
    double Im_x[NMAX];
    double abs_H[NMAX];
    double arg_H[NMAX];
    double f[NMAX];
    float angle;
    double   K, fs,T,finp,Tmax;
    char str;
    char key;

    /*
    float heavi( double t)
    {
    if (t<0 ) return 0;
    if (t>=0) return 1;
    }

    float heavi1( int n)
    {
    if (n<0 ) return 0;
    if (n>=0) return 1;
    }*/

    double getarg(float Re_H,float Im_H)
    {

    return atan2(Im_H,Re_H) * 180 / M_PI;
    }


     int _tmain(int argc, _TCHAR* argv[])
    {


    printf ("\n For input  |H(jw)|, fi(w) from file input '1' ");
    printf ("\n For manual  input |H(jw)|, fi(w) input   '0' ");

    //std::cin>>key;
    scanf( "%s", &key);
    if (key=='0' ) { goto label_02;}
    if (key=='1' ) { goto label_01;}

    label_01:
     printf("\n Open afrcoefs.txt for reading AFR coefficients ");
    if ((fp=fopen("afrcoefs.txt","r"))==NULL)  printf ("\n File opening error");
    fscanf(fp,"%le",&fs);
    fscanf(fp,"%d",&N);
    Tmax=N*1/fs;
    for (k=0 ; k<N;k++)
    {
    // f[k]=k*fs/N;
    //fprintf(fp1,"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
     fscanf (fp,"%d",&n);
     fscanf (fp,"%le",&abs_H[k]);
     fscanf (fp,"%le",&arg_H[k]);
     
      fscanf (fp,"%le",&f[k]);
     
     arg_H[k]=arg_H[k]*M_PI/180;
    }
    fclose(fp);
    goto label_03;

    label_02:

    printf("\n Input fs ,Hz ");
    scanf("%le",&fs);
    printf(" Input  number of points ");
    scanf("%d",&N);
    Tmax=N*1/fs;
    printf("\n Tmax =%lf s", Tmax);
    for (k=0 ; k<N;k++)
    {
    f[k]=k*fs/N;

    //fprintf(fp1,"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
          printf("\n F=%le Hz" ,f[k]);
       printf("\n k=%d;    |H(jw)|[k]=",k );
       scanf ( "%le",&abs_H[k]);
          printf(" k=%d;   arg(H(jw))[k]=",k );
          scanf ( "%le",&arg_H[k]);
    //abs_H[k]=1;
    //arg_H[k]=0;
    arg_H[k]=arg_H[k]*M_PI/180;
    }

    label_03:
     
    for (k=0 ; k<N;k++)
    {
    f[k]=k*fs/N;
    printf("\n n=%d;|H(jw)|= %le ;arg(H(jw))= %le deg;f=%lf Hz",k, abs_H[k] ,  arg_H[k]*180/M_PI , f[k]) ;
     
    }


    printf("\n Input 0 for continue ");
    scanf("%s",&key);
    printf("\n OK. Parsing ...");
     
    // H(jw)=integr(0;inf;h(tau)*exp(-j*w*tau) ; d tau)
    //K=2*3.1415926/fs;
    /*
    DFT
    X[k]=sum(n=0;N-1; x[n]*exp(-2*pi*j*k*n/N) ), k=0,...,N-1
    IDFT
    x[n]=(1/N)*sum(k=0;N-1;X[k]*exp(2*pi*k*n/N) ), n=0,...,N-1
    100*10^-6*500
    */

    label1:;


    //cidft
    printf("\n Open hcoefs1.txt for writing coeffs h[n] ");
    if ((fp1=fopen("hcoefs1.txt","a"))==NULL)  printf ("\n File opening error");
    fprintf(fp1,"\n\n  %le   %d", fs, N );
    for (n=0 ; n<N; n++)
    {
     Re_x[n]=0;
     Im_x[n]=0;
      for(k=0;k<N;k++)
      {    
       angle=2*M_PI*k*n/N;   
            //Uoutp[w]/Uinp[w]   exp(j( fi out[w]-fi inp[w]) )
       Re_H[k]=abs_H[k]*cos(arg_H[k]);
       Im_H[k]=abs_H[k]*sin(arg_H[k]); 
        Re_x[n] = Re_x[n]+  Re_H[k]*cos(angle)-Im_H[k]*sin(angle);
        Im_x[n] = Im_x[n]+  Re_H[k]*sin(angle) + Im_H[k]*cos(angle);
      }
     Re_x[n]=Re_x[n]/N ;
     Im_x[n]=Im_x[n]/N ;
     // n scaled CIDFT  , not used in this, use if sygnal synth.
     //Re_x[n]=Re_x[n]*N;  
     //Im_x[n]=Im_x[n]*N;  
        //printf("\n n=%d  Re(h(n))=%lf  Im(h(n))=%lf ",n, Re_x[n] ,Im_x[n]  ) ;
        printf("\n n=%d   h[n]=%lf   ",n, Re_x[n]    ) ;
        fprintf(fp1, "\n  %d   %le   ",n, Re_x[n]  ) ; 
    }

     fclose (fp1);
     printf( "\n  Re(h(n)=h[n]), use it to load into FIR  "   ) ;
    label3:;
    printf("\n Input 0 for exit ");
    scanf("%s",&key);
    return 0;
    }

    // stdafx.cpp : source file that includes just the standard includes
    // afrtohtcoefs.pch will be the pre-compiled header
    // stdafx.obj will contain the pre-compiled type information

    #include "stdafx.h"

    // TODO: reference any additional headers you need in STDAFX.H
    // and not in this file

    // stdafx.h : include file for standard system include files,
    // or project specific include files that are used frequently, but
    // are changed infrequently
    //

    #pragma once


    #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers

    #include <stdio.h>
    #include <tchar.h>
     #include <iostream>
    #include <stdlib.h>
    #include <cmath>

    // TODO: reference additional headers your program requires here

    afrcoefs.txt
       4.410000e+004  16
      0   1.000000    0.000000    0.000000 
      1   1.000000    -168.750000    2756.250000 
      2   1.000000    22.500000    5512.500000 
      3   1.000000    -146.249999    8268.750000 
      4   0.000000    45.000005    11025.000000 
      5   0.000000    -123.750001    13781.250000 
      6   0.000000    67.500002    16537.500000 
      7   0.000000    78.750002    19293.750000 
      8   0.000000    90.000002    22050.000000 
      9   0.000000    -78.749998    24806.250000 
      10   0.000000    -67.499999    27562.500000 
      11   0.000000    123.750002    30318.750000 
      12   0.000000    135.000007    33075.000000 
      13   1.000000    146.250002    35831.250000 
      14   1.000000    -22.499997    38587.500000 
      15   1.000000    168.750003    41343.750000 
       4.410000e+004  16
      0   1.000001e+000    0.000000e+000    0.000000e+000 
      1   1.000000e+000    -1.687500e+002    2.756250e+003 
      2   9.999998e-001    2.250004e+001    5.512500e+003 
      3   1.000000e+000    -1.462500e+002    8.268750e+003 
      4   3.961544e-007    -2.516211e+001    1.102500e+004 
      5   5.013260e-007    5.027288e+001    1.378125e+004 
      6   3.873424e-007    1.400188e+002    1.653750e+004 
      7   4.839081e-007    1.122935e+002    1.929375e+004 
      8   3.130084e-007    1.795795e+002    2.205000e+004 
      9   4.799148e-007    -1.126101e+002    2.480625e+004 
      10   3.857124e-007    -1.408047e+002    2.756250e+004 
      11   4.938736e-007    -5.018234e+001    3.031875e+004 
      12   3.914724e-007    2.715836e+001    3.307500e+004 
      13   1.000000e+000    1.462500e+002    3.583125e+004 
      14   9.999998e-001    -2.250004e+001    3.858750e+004 
      15   1.000000e+000    1.687500e+002    4.134375e+004 
       4.410000e+004  16
      0   1.000001e+000    0.000000e+000    0.000000e+000 
      1   1.000000e+000    -1.687500e+002    2.756250e+003 
      2   9.999998e-001    2.250004e+001    5.512500e+003 
      3   1.000000e+000    -1.462500e+002    8.268750e+003 
      4   3.949392e-007    -2.565670e+001    1.102500e+004 
      5   4.989970e-007    5.024488e+001    1.378125e+004 
      6   3.867226e-007    1.403127e+002    1.653750e+004 
      7   4.821592e-007    1.124314e+002    1.929375e+004 
      8   3.130000e-007    1.800000e+002    2.205000e+004 
      9   4.821592e-007    -1.124314e+002    2.480625e+004 
      10   3.867226e-007    -1.403127e+002    2.756250e+004 
      11   4.989970e-007    -5.024488e+001    3.031875e+004 
      12   3.949392e-007    2.565670e+001    3.307500e+004 
      13   1.000000e+000    1.462500e+002    3.583125e+004 
      14   9.999998e-001    -2.250004e+001    3.858750e+004 
      15   1.000000e+000    1.687500e+002    4.134375e+004 

    hcoefs1.txt

      4.410000e+004   16
      0   -4.854692e-002  
      1   3.078799e-002  
      2   6.781642e-002  
      3   -7.925028e-003  
      4   -9.804488e-002  
      5   -3.848729e-002  
      6   1.898830e-001  
      7   4.045169e-001  
      8   4.045168e-001  
      9   1.898829e-001  
      10   -3.848736e-002  
      11   -9.804485e-002  
      12   -7.924644e-003  
      13   6.781681e-002  
      14   3.078794e-002  
      15   -4.854676e-002  

    • Edited by USERPC01 Tuesday, January 23, 2018 10:03 PM adding afrcoefs.txt
    Tuesday, January 23, 2018 9:58 PM

All replies

  • Program for converting h(t) coefficients of the FIR  to AFR  data

    // irtoafr.cpp : Defines the entry point for the console application.
    //

    #include "stdafx.h"

    #include <cstdio> // or  cstdio.h
    #include <cstdlib> //or  cstdlib.h
    #include <cmath>   //or  cmath.h

    #define M_PI 3.1415926535897932384626433832795
    #define PI M_PI // 3.14159265
    #define  NMAX 512
    FILE *fp, *fp1;
    int w,t,N;
    double h[NMAX];
    double h_im[NMAX];
    double Re_H[NMAX];
    double Im_H[NMAX];
    double Re_x[NMAX];
    double Im_x[NMAX];
    double abs_H[NMAX];
    double arg_h[NMAX];
    double f[NMAX];
    double angle;
    double   K, fs,T,Tmax;
    char str;
    char key;
    int k,n;


    /*
    float x[N];
    float y[N];

    double my_FIR(double sample_data, int Nf)
    {
      double result = 0;
      for ( int i = Nf - 2 ; i >= 0 ; i-- )
      {
        x[i + 1] = x[i];
        y[i + 1] = y[i];
      }
      x[0] = (double)sample_data;
      for (int k = 0; k < N; k++)
      {
        result = result + x[k]*h[k];
      }
      y[0] = result;
      return ((double)result);
    }

    */


    float heavi( double t)
    {
    if (t<0 ) return 0;
    if (t>=0) return 1;
    }

    float heavi1( int n)
    {
    if (n<0 ) return 0;
    if (n>=0) return 1;
    }


    float dirac( int n)
    {
    if (n!=0 ) return 0;
    if (n=0) return 1;
    }

    double getarg(float Re_H,float Im_H)
    {
    return atan2(Im_H,Re_H) * 180/ PI;
    }


    int _tmain(int argc, _TCHAR* argv[])
    {
    //double tau=0.02;
    //printf(" Input tau , s ");
    //scanf("%le",&tau);
    printf ("\n For input  h(t) from file input '1' ");
    printf ("\n For manual  input  h(t) input   '0' ");
    scanf("%s",&key);
    if (key=='0' ) { goto label_02;}
    if (key=='1' ) { goto label_01;}

    label_01:
      printf("\n Open htkoefs.txt for reading h[n] coefficients ");
    if ((fp=fopen("htkoefs.txt","r"))==NULL)  printf ("\n File opening error");
    fscanf(fp,"%le",&fs);
    fscanf(fp,"%d",&N);
    Tmax=N*1/fs;
    for (n=0 ; n<N;n++)
    {
     fscanf (fp,"%d",&k);
     fscanf (fp,"%le",&h[n]);
    }
    fclose(fp);
    goto label_03; 
    label_02:
     
     
    printf("\n Input fs ,Hz ");
    scanf("%le",&fs);
    printf(" Input  number of points ");
    scanf("%d",&N);
    Tmax=N*1/fs;
    printf("\n Tmax =%lf s \n", Tmax); 
     
    for (n=0 ; n<N;n++)
    {
     T=n*1/(N*fs); 
     printf("\nT=%le s; n =%d ; input h(n): ",T,n);
     scanf ("%le",&h[n]);

    /*
    h[0]=-0.0485469;
    h[1]=+0.0307880;
    h[2]=+0.0678165;
    h[3]=-0.0079250;
    h[4]=-0.0980449;
    h[5]=-0.0384873;
    h[6]=+0.1898828;
    h[7]=+0.4045168;
    h[8]=+0.4045168;
    h[9]=+0.1898828;
    h[10]=-0.0384873;
    h[11]=-0.0980449;
    h[12]=-0.0079250;
    h[13]=+0.0678165;
    h[14]=+0.0307880;
    h[15]=-0.0485469;
    */
    //tau=0.1*Tmax;
    //h[t]=1-exp(-T/tau);
    }

    label_03:
    printf("\n Fs=%lf Hz",fs);
    printf("\n N=%d  points",N);
    for (n=0 ; n<N;n++)
    {
     T=n*1/(N*fs); 
     printf("\n t=%le s; h[t]=%le   ",T,h[n]);
    }

    printf("\n Input 0 for next ");
    scanf("%s",&key);
    printf("\n OK. Parsing ... ");
     
    // H(jw)=integr(0;inf;h(tau)*exp(-j*w*tau) ; d tau)
    //
    /*
    DFT
    X[k]=sum(n=0;N-1; x[n]*exp(-2*pi*j*k*n/N) ), k=0,...,N-1
    IDFT
    x[n]=(1/N)*sum(k=0;N-1;X[k]*exp(2*pi*k*n/N) ), n=0,...,N-1
    */

    /*********************************************************/

    label1:
     printf("\n Open afrcoefs.txt for reading AFR coefficients ");
    if ((fp1=fopen("afrcoefs.txt","a"))==NULL)  printf ("\n File opening error");
    fprintf(fp1,"\n   %le  %d ",fs,N);
    /*
    fprintf(fp1,"\n Re(H(jw))  ; Im(H(jw)) ");
    for (k=0 ; k<N; k++)
    {
     Re_H[k]=0;
     Im_H[k]=0;
      for(n=0;n<N;n++)
      {    
       angle=-2*PI*k*n/N;
     //  Re_H[k] = Re_H[k]+ h[n]*cos(angle) - h_im[n]*sin(angle);
     //  Im_H[k] = Im_H[k]+ h[n]*sin(angle) + h_im[n]*cos(angle);
         Re_H[k] = Re_H[k]+ h[n]*cos(angle) ;//- h_im[n]*sin(angle);
         Im_H[k] = Im_H[k]+ h[n]*sin(angle) ;//+ h_im[n]*cos(angle);
      }

    // 1/n scaled DFT(h(t)*1(t))
     //Re_H[k]=Re_H[k]/N ;
     //Im_H[k]=Im_H[k]/N ;
       abs_H[k]=pow(((Re_H[k]*Re_H[k])+(Im_H[k]*Im_H[k])), 0.5);
       arg_h[k]=getarg(Re_H[k],Im_H[k]);
       f[k]=k*fs/N; 
         printf("\n n=%d Re(H(jw))= %lf ; Im(H(jw))= %lf ; f=%lf Hz",k, Re_H[k] ,Im_H[k], f[k]) ;
         fprintf(fp1,"\n n=%d Re(H(jw))= %lf ; Im(H(jw))= %lf ; f=%lf Hz",k, Re_H[k] ,Im_H[k], f[k]) ;
    }

    printf("\n Input 0 for continue ");
    scanf("%s",&key);
    */
     //fprintf(fp1,"\n |H(jw)| , arg( H(jw) ) ");
     
    for (k=0 ; k<N; k++)
    {
     Re_H[k]=0;
     Im_H[k]=0;
      for(n=0;n<N;n++)
      {   
       angle=-2*PI*k*n/N;
     //  Re_H[k] = Re_H[k]+ h[n]*cos(angle) - h_im[n]*sin(angle);
     //  Im_H[k] = Im_H[k]+ h[n]*sin(angle) + h_im[n]*cos(angle);
         Re_H[k] = Re_H[k]+ h[n]*cos(angle) ;//- h_im[n]*sin(angle);
         Im_H[k] = Im_H[k]+ h[n]*sin(angle) ;//+ h_im[n]*cos(angle);
      }

     // 1/n scaled DFT(h(t)*1(t))
     //Re_H[k]=Re_H[k]/N ;
     //Im_H[k]=Im_H[k]/N ;
     
       abs_H[k]=pow(((Re_H[k]*Re_H[k])+(Im_H[k]*Im_H[k])), 0.5);
       arg_h[k]=getarg(Re_H[k],Im_H[k]);
       f[k]=k*fs/N;
       printf("\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
        //fprintf(fp1,"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
       fprintf(fp1,"\n  %d   %le    %le    %le  ",k, abs_H[k] ,arg_h[k], f[k]) ;
    }
    fclose(fp1);
    printf("\n Input 0 for exit ");
    scanf("%s",&key);

    return 0;
    }

    //cidft
    /*
    for (n=0 ; n<N; n++)
    {
     Re_x[n]=0;
     Im_x[n]=0;
      for(k=0;k<N;n++)
      {  
    //  angle = 2 * PI * k * N / numidft
    //  idft_xnre(N + 1) = idft_xnre(N + 1) + idft_Xkre(k + 1) * Cos(angle) - idft_Xkim(k + 1) * Sin(angle)
    //  idft_xnim(N + 1) = idft_xnim(N + 1) + idft_Xkre(k + 1) * Sin(angle) + idft_Xkim(k + 1) * Cos(angle) 
       angle=2*PI*k*n/N;
       Re_H[k]=abs_H[k]*cos(angle);
       Im_H[k]=abs_H[k]*sin(angle);
        Re_x[n] = Re_x[n]+ Re_H[k]*cos(angle)- Im_H[k]*sin(angle);
        Im_x[n] = Im_x[n]+ Re_H[k]*sin(angle)+ Im_H[k]*cos(angle);
      }
     Re_x[n]=Re_x[n]/N ;
     Im_x[n]=Im_x[n]/N ;
     // n scaled CIDFT( )
     //Re_x[n]=Re_x[n]*N;  
     //Im_x[n]=Im_x[n]*N;  
     //  abs_H[k]=pow(((Re_H[k]*Re_H[k])+(Im_H[k]*Im_H[k])), 0.5);
     //  arg_h[k]=getarg(Re_H[k],Im_H[k]);
     //  T[k]=k*fs/N;
       printf("\n n=%d  Re(h(n))=%d   ",n, Re_x[n] ,Im_x[n]  ) ;

    }

    */

    // stdafx.cpp : source file that includes just the standard includes
    // irtoafr.pch will be the pre-compiled header
    // stdafx.obj will contain the pre-compiled type information

    #include "stdafx.h"

    // TODO: reference any additional headers you need in STDAFX.H
    // and not in this file

    // stdafx.h : include file for standard system include files,
    // or project specific include files that are used frequently, but
    // are changed infrequently
    //

    #pragma once


    #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
    #include <stdio.h>
    #include <tchar.h>

    // TODO: reference additional headers your program requires here

    Tuesday, January 23, 2018 10:01 PM
  • Implementation of the FIR

    // fir.cpp : Defines the entry point for the console application.
    //

    #include "stdafx.h"

    #include <cstdio> // or  cstdio.h
    #include <cstdlib> //or  cstdlib.h
    #include <cmath>   //or  cmath.h

     
    #define  NMAX 512
    FILE *fp, *fp1, *fp2;
    int w,t,N,N1;
    double h[NMAX];
    double xin[NMAX];
    double yout[NMAX];
    double f[NMAX];
     
    double   K, fs,T,Tmax;
    char str;
    char key;
    int k,n,i;


    double my_FIR(double sample_data, int Nf,  double *h)
    {
    static double x[NMAX];
    static double y[NMAX];

      double result = 0;
      for ( int i = Nf - 2 ; i >= 0 ; i-- )
      {
        x[i + 1] = x[i];
        y[i + 1] = y[i];
      }
      x[0] = (double)sample_data;
      for (int k = 0; k < Nf; k++)
      {
        result = result + x[k]*h[k];
      }
      y[0] = result;
      return ((double)result);
    }

     


    float heavi( double t)
    {
    if (t<0 ) return 0;
    if (t>=0) return 1;
    }

    float heavi1( int n)
    {
    if (n<0 ) return 0;
    if (n>=0) return 1;
    }


    float dirac( int n)
    {
    if (n!=0 ) return 0;
    if (n=0) return 1;
    }

     


    int _tmain(int argc, _TCHAR* argv[])
    {
     
    printf ("\n For input  h(t) from file input '1' ");
    printf ("\n For manual  input  h(t) input   '0' ");
    scanf("%s",&key);
    if (key=='0' ) { goto label_02;}
    if (key=='1' ) { goto label_01;}

    label_01:

     printf("\n Open htkoefs.txt for reading h[n] coefficients ");
    if ((fp=fopen("htkoefs.txt","r"))==NULL)  printf ("\n File opening error");
    fscanf(fp,"%le",&fs);
    fscanf(fp,"%d",&N);

    for (n=0 ; n<N;n++)
    {
     fscanf (fp,"%d",&k);
     fscanf (fp,"%le",&h[n]);
    }


    fclose(fp);
    goto label_03;
     
    label_02:
     
     
    printf("\n Input fs ,Hz " );
    scanf("%le",&fs);
    printf(" Input order ");
    scanf("%d",&N);
      
     
    for (n=0 ; n<N;n++)
    {
     
     printf("\n   n =%d ; input h[n]: " ,n);
     scanf ("%le",&h[n]);

    }

    label_03:
    printf("\n Fs coefs=%lf Hz ",fs);
    printf("\n order=  %d     ",N);

    for (n=0 ; n<N;n++)
    {
     printf("\n  n =%d ; h[n]=%le   ",n,h[n]);
    }

    printf("\n Input 0 for continue ");
    scanf("%s",&key);
    printf("\n OK.   ");


    printf ("\n For input  x(t) from file input '1' ");
    printf ("\n For manual  input  x(t) input   '0' ");
    scanf("%s",&key);
    if (key=='1' ) { goto label_05;}
    if (key=='0' ) { goto label_06;}

    label_05:

     printf("\n Open input.txt for reading x[n]   ");
    if ((fp1=fopen("input.txt","r"))==NULL)  printf ("\n File opening error ");
    fscanf(fp1,"%le",&fs);
    fscanf(fp1,"%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     fscanf (fp1,"%d",&k);
     fscanf (fp1,"%le",&xin[n]);
    }


    fclose(fp1);

    goto label_07;

    label_06:

    printf("\n Input fs,Hz , ( %lf ) ",fs) ;
    scanf( "%le",&fs);
    printf("\n Input number of samples of xin[n]  ") ;
    scanf( "%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     T=n*1/(N1*fs); 
     printf ("\n n=%d ; t=%le s; xin[n]= ",n,T);
     scanf ( "%le",&xin[n]);
    }

     label_07:
    printf( "\n fs=%le Hz ",fs);
    printf( "\n Nsamples= %d ",N1);

    for (n=0 ; n<N1;n++)
    {
     T=n*1/(N1*fs); 
     printf("\n n=%d x[n]=%le,  t=%le s",n,xin[n],T);
    }

    printf("\n Input 0 for continue ");
    scanf("%s",&key);
    printf("\n OK.   ");

    /*********************************************************/

    label1:
      printf("\n Open output.txt for writing sequence   ");

    if ((fp2=fopen("output.txt","a"))==NULL)  printf ("\n File opening error");
    fprintf(fp2,"\n   %le  %d \n",fs,N);

    for (i=0 ;  i<N1;i++)
    {
     T=i*1/(N1*fs);
    yout[i]= my_FIR(xin[i], N,h);
     printf( "\n n=%d t=%le s; yout[n]= %le ",i,T,yout[i] );
    fprintf(fp2,"\n %d  %le ",i,yout[i] );
    }


    fclose(fp2);
    printf("\n Input 0 for exit ");
    scanf("%s",&key);

    return 0;
    }

     

    // stdafx.cpp : source file that includes just the standard includes
    // fir.pch will be the pre-compiled header
    // stdafx.obj will contain the pre-compiled type information

    #include "stdafx.h"

    // TODO: reference any additional headers you need in STDAFX.H
    // and not in this file

    // stdafx.h : include file for standard system include files,
    // or project specific include files that are used frequently, but
    // are changed infrequently
    //

    #pragma once


    #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
    #include <stdio.h>
    #include <tchar.h>

    // TODO: reference additional headers your program requires here

    input.txt
       4.410000e+004   16

     0  1.000000e+001  
     1  7.071068e+000  
     2  6.123032e-016  
     3  -7.071068e+000  
     4  -1.000000e+001  
     5  -7.071068e+000  
     6  -1.836910e-015  
     7  7.071068e+000  
     8  1.000000e+001  
     9  7.071068e+000  
     10  3.061516e-015  
     11  -7.071068e+000  
     12  -1.000000e+001  
     13  -7.071068e+000  
     14  -4.286122e-015  
     15  7.071068e+000 

    htkoefs.txt

      4.410000e+004   16
      0   0 
      1   1 
      2   0  
      3   0  
      4   0  
      5   0  
      6   0 
      7   0  
      8   0  
      9   0  
      10  0  
      11  0  
      12  0  
      13  0  
      14  0  
      15  0 

    output.txt


       4.410000e+004  16

     0  0.000000e+000
     1  1.000000e+001
     2  7.071068e+000
     3  6.123032e-016
     4  -7.071068e+000
     5  -1.000000e+001
     6  -7.071068e+000
     7  -1.836910e-015
     8  7.071068e+000
     9  1.000000e+001
     10  7.071068e+000
     11  3.061516e-015
     12  -7.071068e+000
     13  -1.000000e+001
     14  -7.071068e+000
     15  -4.286122e-015


    • Edited by USERPC01 Tuesday, January 23, 2018 10:08 PM fixed bugs
    Tuesday, January 23, 2018 10:07 PM
  • Program for spectrum calculating, using CDFT,CIDFT, Mathlab's algorithm for Hilbert transform

    You can use complex.h, vector.h for data arrays and pointered arrays for arrays.

    How to use radix more then 20000 in  #define NMAX  , how to fix  memory mapping error in arrays if  NMAX>20000...32768?

    // getspec.cpp : Defines the entry point for the console application.
    //

    #include "stdafx.h"

    #include <cstdio>
    #include <cstdlib>
    #include <cmath>
    #include <limits>
    #include <iostream>
    using namespace std;

    #define NMAX 2048

    #define M_PI 3.1415926535897932384626433832795

    typedef struct {

    double re[NMAX];
    double im[NMAX];
    } CData;

     typedef struct  {
         double data[NMAX];
     } Array;


    CData  cdft(CData xn  ,int N )
     {
       
     
     CData Xk ;
     int k,n  ;
     double angle;

     for (k=0;k<N;k++)
      {
      Xk.re[k]=0;
      Xk.im[k]=0;
      for(n=0;n<N;n++)
        {
        angle=-2*M_PI*k*n/N;  
        Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
        Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
        }

     }

     return Xk;
     }

     CData  cidft(CData idft_Xk, int N  )
     {

     int k,n ;
     double angle;
     CData  idft_xn;
     
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (n=0;n<N;n++)
      {
      idft_xn.re[n]=0;
      idft_xn.im[n]=0;
      for(k=0;k<N;k++)
       {
     
       angle=2*M_PI*k*n/N;
       idft_xn.re[n]=idft_xn.re[n]+idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
       idft_xn.im[n]=idft_xn.im[n]+idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
       }

     idft_xn.im[n]=idft_xn.im[n]/N;
     idft_xn.re[n]=idft_xn.re[n]/N;

     }

    return idft_xn ;
     }

    CData ScaleDFT(CData xk,int  num)
     {

     int i;
     for (i=0; i<num ; i++)
       {
        xk.re[i]=xk.re[i]/num;
        xk.im[i]=xk.im[i]/num;
       }
     return xk;
     }

    CData ScaleIQidft(CData Xk, int num)
     {
     int i;
       for (i=0;  i<num ; i++)
      {
       Xk.re[i]=Xk.re[i]*num;
       Xk.im[i]=Xk.im[i]*num;
      }

    return Xk;

    }

     CData hilbert(double *xr , int N)
     {
      int  k,n,i, no2 ;
      double angle;
      double h[NMAX];
      double dft__xn_re[NMAX];
      double dft__xn_im[NMAX];
       double dft__Xk_re[NMAX];
       double dft__Xk_im[NMAX];
         double idft__XK_re[NMAX];
         double idft__XK_im[NMAX];
       CData idft__XN;

     
      if (N==0)
      {
      printf("\n Error");
      idft__XN.re[0]=0;
      idft__XN.im[0]=0;
      return idft__XN;
      }
      //xnreal=real(xr);
      no2=(int) floor ((double)N/2);
     
      for (i=0;i<N;i++)
      {
       dft__xn_re[i]=xr[i];
       dft__xn_im[i]=0;
      }

     //fft
     // dft__Xk=dft(dft__xn,n);

    for (k=0;k<N;k++)
      {
      dft__Xk_re[k]=0;
      dft__Xk_im[k]=0;
      for(n=0;n<N;n++)
        {
        angle=-2*M_PI*k*n/N;
        dft__Xk_re[k]+= dft__xn_re[n]*cos(angle)- dft__xn_im[n]*sin(angle);
        dft__Xk_im[k]+= dft__xn_re[n]*sin(angle)+ dft__xn_im[n]*cos(angle);
        }

     }

     //X={Xk.re,Xk.im};

     for (i=0;i<=N;i++)
      {
      h[i]=0;
      }


      if((N>0)&&(( 2*floor((double)N/2) )==(N))) // n is even
      {

      h[0]=1;
       h[N/2]=1;
       for(k=1;k<N/2;k++)
       {
       h[k]=2;
       }
      }
      else //n is odd
      {
        if(N>0)
        {
         h[0]=1;
         for(k=1;k<(N+1)/2;k++)
         {
         h[k]=2;
         }
        }
      }

     for(i=0;i<N;i++)
      {
      idft__XK_re[i]= dft__Xk_re[i]*h[i];
      idft__XK_im[i]= dft__Xk_im[i]*h[i];
      }


     //idft__XN= idft(idft__XK, n);
     
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (n=0;n<N;n++)
      {
      idft__XN.re[n]=0;
      idft__XN.im[n]=0;
      for(k=0;k<N;k++)
       {
       //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
       angle=2*M_PI*k*n/N;
       idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
       idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
       }
     idft__XN.im[n]=idft__XN.im[n]/N;
     idft__XN.re[n]=idft__XN.re[n]/N;

     }


    return idft__XN;
     }

    CData getanalyticsignal(double *xn, int num)
     {
     CData x_n,Hilb;
     int i;
       //x=x+jy,y=0
       //hilb(x)=a+jb
       //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
       //I=x(t), Q=j*Hilb(xre(t))

    Hilb=hilbert(xn,num);
     
     for (i=0; i<num; i++)
     {
      x_n.re[i]=xn[i];
      x_n.im[i]=Hilb.im[i];
     }
     
     return x_n;
     }


     
    double Arg(double re, double im)
    {
     
    return atan2 (im,re)  ; 
    }

    CData  get_mod_fi(CData x, int n, int mode   )
     {
      int i;
      CData res;
      double fi;
      for (i=0;i<n;i++)
      {
      res.re[i]=pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
      //double treshh=1e9;
      // x.im[i]=x.im[i]*treshh;
      // x.im[i]=floor(x.im[i])/treshh;
      // x.re[i]=x.re[i]*treshh;
      // x.re[i]=floor(x.re[i])/treshh;
      res.im[i]=Arg(x.re[i],x.im[i]);
      if (mode==1)  { res.im[i]=res.im[i]*180/M_PI;  }
              
      }
      return  res;
     }

    CData GetSpectrum(double *xn, int num)
    {
     
    CData xn_iq, Xf_IQ, Xf_modfi;
    xn_iq=getanalyticsignal( xn,  num);
     
    Xf_IQ=cdft(xn_iq  ,num );
    Xf_IQ=ScaleDFT(Xf_IQ,num);
    Xf_modfi=get_mod_fi(Xf_IQ, num,1   );
     
    return Xf_modfi;
    }

    Array getenvelope(double *x, int num)
     {
     int i;
     CData  y;
     Array env;
     y=hilbert(x,num);
     for (i=0;i<num;i++)
      {
      env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
      }
      return env;
     }

    FILE   *fp1,*fp2 ;


    int _tmain(int argc, _TCHAR* argv[])
    {

    double xin[NMAX];
    CData  Xf;

    char key;
    int k,n,i,N1;
    double fs,T,x,f;


    printf ("\n For input  x(t) from file input '1' ");
    printf ("\n For manual  input  x(t) input   '0' ");
    scanf("%s",&key);
    if (key=='1' ) { goto label_05;}
    if (key=='0' ) { goto label_06;}

    label_05:

     printf("\n Open xinput.txt for reading x[n]   ");
    if ((fp1=fopen("xinput.txt","r"))==NULL)  printf ("\n File opening error ");
    fscanf(fp1,"%le",&fs);
    fscanf(fp1,"%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     fscanf (fp1,"%d",&k);
     fscanf (fp1,"%le",&xin[n]);
    }


    fclose(fp1);

    goto label_07;

    label_06:

    printf("\n Input fs,Hz , ( %lf ) ",fs) ;
    scanf( "%le",&fs);
    printf("\n Input number of samples of xin[n]  ") ;
    scanf( "%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     T=n*1/(N1*fs); 
     printf ("\n n=%d ; t=%le s; xin[n]= ",n,T);
     scanf ( "%le",&xin[n]);
    }

     label_07:
    printf( "\n fs=%le Hz ",fs);
    printf( "\n Nsamples= %d ",N1);

    for (n=0 ; n<N1;n++)
    {
     T=n*1/(N1*fs); 
     printf("\n n=%d x[n]=%le,  t=%le s",n,xin[n],T);
    }

    printf("\n Input 0 for continue ");
    scanf("%s",&key);
    printf("\n OK.   ");


    /*****************************/

    label1:
    //Xf.re= |Xf|,  Xf.im -arg(Xf) (only in this program)
    printf("\n Obtaining spectrum of x[n]");
    Xf=GetSpectrum(xin, N1);

    printf("\n Spectrum of x[n] : ");

     if ((fp2=fopen("spectrum.txt","a"))==NULL)  printf ("\n File opening error");
     fprintf(fp2,"\n   %le   %d \n",fs,N1);
     //fprintf(fp2, "\n n;\t f, Hz;\t |Xf(jw)| ;\t  arg(Xf), deg ");
    for (i=0; i<N1;i++)
    {
    f=fs*i/N1;
    //cout <<"\n f="<< f<<" Hz;|Xf(jw)|="<<Xf.re[i]<<";arg(Xf)="<<Xf.im[i]<<" deg ";
     printf ("\n f=%le Hz;|Xf(jw)|=%le;arg(Xf)=%lf deg ",  f,Xf.re[i], Xf.im[i]);
     fprintf(fp2,"\n %d   %le  %le",i ,Xf.re[i],Xf.im[i]);

    }
    printf ("\n  input   '0' for exit ");
    char key1;
    cin>>key1;

    return 0;

    }

    // stdafx.cpp : source file that includes just the standard includes
    // getspec.pch will be the pre-compiled header
    // stdafx.obj will contain the pre-compiled type information

    #include "stdafx.h"

    // TODO: reference any additional headers you need in STDAFX.H
    // and not in this file

    // stdafx.h : include file for standard system include files,
    // or project specific include files that are used frequently, but
    // are changed infrequently
    //

    #pragma once


    #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
    #include <stdio.h> //edit
    #include <tchar.h> //edit

    // TODO: reference additional headers your program requires here

    xinput.txt

       4.410000e+004   16

     0  1.000000e+001  
     1  7.071068e+000  
     2  6.123032e-016  
     3  -7.071068e+000  
     4  -1.000000e+001  
     5  -7.071068e+000  
     6  -1.836910e-015  
     7  7.071068e+000  
     8  1.000000e+001  
     9  7.071068e+000  
     10  3.061516e-015  
     11  -7.071068e+000  
     12  -1.000000e+001  
     13  -7.071068e+000  
     14  -4.286122e-015  
     15  7.071068e+000   

    spectrum.txt

       4.410000e+004   16

     0   7.216450e-016  -1.126199e+002
     1   7.237769e-016  -8.560129e+001
     2   1.000000e+001  9.223608e-015
     3   9.436896e-016  1.800000e+002
     4   3.373574e-016  3.464047e+001
     5   5.950643e-015  -5.660642e+001
     6   1.330312e-007  1.800000e+002
     7   7.196406e-015  -2.652734e+000
     8   3.760637e-015  -7.363206e+001
     9   3.624946e-015  -1.272209e+002
     10   2.820163e-015  1.438068e+002
     11   5.525033e-015  -1.015922e+002
     12   4.814310e-015  7.423349e+000
     13   1.110917e-014  6.612472e+001
     14   6.445273e-015  -9.246812e+001
     15   8.760391e-015  -1.326889e+002



    • Edited by USERPC01 Tuesday, January 23, 2018 10:29 PM adding xinput.txt
    Tuesday, January 23, 2018 10:14 PM
  • Program for signal synth. from freq. data

    // synthsig.cpp : Defines the entry point for the console application.
    //

    #include "stdafx.h"
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <iostream>

    #define NMAX 1000
    #define M_PI 3.1415926535897932384626433832795

    typedef struct {
    double re[NMAX];
    double im[NMAX];
    } CData;

     typedef struct  {
         double data[NMAX];
     } Array;


    CData  cdft(CData xn  ,int N )
     {
       
      
     CData Xk ;
     int k,n  ;
     double angle;

     for (k=0;k<N;k++)
      {
      Xk.re[k]=0;
      Xk.im[k]=0;
      for(n=0;n<N;n++)
        {
        angle=-2*M_PI*k*n/N;  
        Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
        Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
        }

     }

     return Xk;
     }

     CData  cidft(CData idft_Xk, int N  )
     {

     int k,n ;
     double angle;
     CData  idft_xn;
     
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (n=0;n<N;n++)
      {
      idft_xn.re[n]=0;
      idft_xn.im[n]=0;
      for(k=0;k<N;k++)
       {
     
       angle=2*M_PI*k*n/N;
       idft_xn.re[n]=idft_xn.re[n]+idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
       idft_xn.im[n]=idft_xn.im[n]+idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
       }

     idft_xn.im[n]=idft_xn.im[n]/N;
     idft_xn.re[n]=idft_xn.re[n]/N;

     }

    return idft_xn ;
     }

    CData ScaleDFT(CData xk,int  num)
     {

     int i;
     for (i=0; i<num ; i++)
       {
        xk.re[i]=xk.re[i]/num;
        xk.im[i]=xk.im[i]/num;
       }
     return xk;
     }

    CData ScaleIQidft(CData Xk, int num)
     {
     int i;
       for (i=0;  i<num ; i++)
      {
       Xk.re[i]=Xk.re[i]*num;
       Xk.im[i]=Xk.im[i]*num;
      }

    return Xk;

    }

     CData hilbert(double *xr , int N)
     {
      int  k,n,i, no2 ;
      double angle;
      double h[NMAX];
      double dft__xn_re[NMAX];
      double dft__xn_im[NMAX];
       double dft__Xk_re[NMAX];
       double dft__Xk_im[NMAX];
         double idft__XK_re[NMAX];
         double idft__XK_im[NMAX];
       CData idft__XN;

     
      if (N==0)
      {
      printf("\n Error");
      idft__XN.re[0]=0;
      idft__XN.im[0]=0;
      return idft__XN;
      }
      //xnreal=real(xr);
      no2=(int)floor ((double)N/2);
     
      for (i=0;i<N;i++)
      {
       dft__xn_re[i]=xr[i];
       dft__xn_im[i]=0;
      }

     //fft
     // dft__Xk=dft(dft__xn,n);

    for (k=0;k<N;k++)
      {
      dft__Xk_re[k]=0;
      dft__Xk_im[k]=0;
      for(n=0;n<N;n++)
        {
        angle=-2*M_PI*k*n/N;
        dft__Xk_re[k]+= dft__xn_re[n]*cos(angle)- dft__xn_im[n]*sin(angle);
        dft__Xk_im[k]+= dft__xn_re[n]*sin(angle)+ dft__xn_im[n]*cos(angle);
        }

     }

     //X={Xk.re,Xk.im};

     for (i=0;i<=N;i++)
      {
      h[i]=0;
      }


      if((N>0)&&((2*floor((double) N/2))==(N))) // n is even
      {

      h[0]=1;
       h[N/2]=1;
       for(k=1;k<N/2;k++)
       {
       h[k]=2;
       }
      }
      else //n is odd
      {
        if(N>0)
        {
         h[0]=1;
         for(k=1;k<(N+1)/2;k++)
         {
         h[k]=2;
         }
        }
      }

     for(i=0;i<N;i++)
      {
      idft__XK_re[i]= dft__Xk_re[i]*h[i];
      idft__XK_im[i]= dft__Xk_im[i]*h[i];
      }


     //idft__XN= idft(idft__XK, n);
     
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (n=0;n<N;n++)
      {
      idft__XN.re[n]=0;
      idft__XN.im[n]=0;
      for(k=0;k<N;k++)
       {
       //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
       angle=2*M_PI*k*n/N;
       idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
       idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
       }
     idft__XN.im[n]=idft__XN.im[n]/N;
     idft__XN.re[n]=idft__XN.re[n]/N;

     }


    return idft__XN;
     }

    CData getanalyticsignal(double *xn, int num)
     {
     CData x_n,Hilb;
     int i;
       //x=x+jy,y=0
       //hilb(x)=a+jb
       //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
       //I=x(t), Q=j*Hilb(xre(t))

    Hilb=hilbert(xn,num);
     
     for (i=0; i<num; i++)
     {
      x_n.re[i]=xn[i];
      x_n.im[i]=Hilb.im[i];
     }
     
     return x_n;
     }


     CData GetIQfromR_fi(CData Xinp, int numb)
     {
     int i;
     CData res;
      for (i=0;i<numb;i++)
      {
       Xinp.im[i]=Xinp.im[i]*M_PI/180;
      res.re[i]= Xinp.re[i]*cos(Xinp.im[i]);
      res.im[i]= Xinp.re[i]*sin(Xinp.im[i]);
      }

     return res;
     }

    double Arg(double re, double im)
    {
     
    return atan2(im,re); 
    }

    CData  get_mod_fi(CData x, int n, int mode   )
     {
      int i;
      CData res;
      double fi;
      for (i=0;i<n;i++)
      {
      res.re[i]=pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
      //double treshh=1e9;
      // x.im[i]=x.im[i]*treshh;
      // x.im[i]=floor(x.im[i])/treshh;
      // x.re[i]=x.re[i]*treshh;
      // x.re[i]=floor(x.re[i])/treshh;
      res.im[i]=Arg(x.re[i],x.im[i]);
      if (mode==1)  { res.im[i]=res.im[i]*180/M_PI;  }
              
      }
      return  res;
     }

    CData GetSpectrum(double *xn, int num)
    {
     
    CData xn_iq, Xf_IQ, Xf_modfi;
    xn_iq=getanalyticsignal( xn,  num);
     
    Xf_IQ=cdft(xn_iq  ,num );
    Xf_IQ=ScaleDFT(Xf_IQ,num);
    Xf_modfi=get_mod_fi(Xf_IQ, num,1   );
     
    return Xf_modfi;
    }

    CData SynthSignal(CData Xf_modfi, int num)

    {
    CData Xf_IQ,xn_iq;
    Xf_IQ=GetIQfromR_fi( Xf_modfi, num);
    xn_iq=cidft(Xf_IQ,num  );
    xn_iq=ScaleIQidft(xn_iq,  num);
    return xn_iq;
    }

    Array getenvelope(double *x, int num)
     {
     int i;
     CData  y;
     Array env;
     y=hilbert(x,num);
     for (i=0;i<num;i++)
      {
      env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
      }
      return env;
     }

    FILE   *fp1,*fp2 ;


    //int main ()
    int _tmain(int argc, _TCHAR* argv[])
    {

    double xin[NMAX];
    CData  Xf;
    CData xt;
    char key;
    int k,n,i,N1;
    double fs=44100, T,x,f;


    printf ("\n For input  x(t) from file input '1' ");
    printf ("\n For manual  input  x(t) input   '0' ");
    scanf("%s",&key);
    if (key=='1' ) { goto label_05;}
    if (key=='0' ) { goto label_06;}

    label_05:

     printf("\n Open spectrum.txt for reading Xf   ");
    if ((fp1=fopen("spectrum.txt","r"))==NULL)  printf ("\n File opening error ");
    fscanf(fp1,"%le",&fs);
    fscanf(fp1,"%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     fscanf (fp1,"%d",&k);
     fscanf (fp1,"%le", &Xf.re[n]);
     fscanf (fp1,"%le", &Xf.im[n]);
    }


    fclose(fp1);

    goto label_07;

    label_06:

    printf("\n Input fs,Hz , ( %lf ) ", fs) ;
    scanf( "%le",&fs);
    printf("\n Input number of samples of xin[n]  ") ;
    scanf( "%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     
     f=fs*n/N1; 
     printf ("\n n=%d ; f=%le Hz; ",n,f);
     printf ("\n Input |X(jw)| :" );
     scanf ( "%le",&Xf.re[n]);
     printf ("\nInput arg(X(jw)), deg :" );
     scanf ( "%le",&Xf.im[n]);
    }

     label_07:
    printf( "\n fs=%le Hz ",fs);
    printf( "\n Nsamples= %d ",N1);

    for (n=0 ; n<N1;n++)
    {
      f=fs*n/N1;
     printf ("\n f=%le Hz;|Xf(jw)|=%le;arg(Xf)=%le deg", f,Xf.re[n],Xf.im[n]);
    }

    printf("\n Input 0 for continue ");
    scanf("%s",&key);
    printf("\n OK.   ");


    /*****************************/

    label1:
     
    printf("\n Parsing   Xf->x[n]");
    xt=SynthSignal( Xf , N1);

    printf("\n Signal (xinput.txt): ");
     
     
    if ((fp2=fopen("xinput.txt","a"))==NULL)  printf ("\n File opening error");
    fprintf(fp2,"\n   %le   %d \n",fs,N1);

    for (i=0; i<N1;i++)
    {
     
    printf ("\n  %d ;   %le ;   %le ;",  i,xt.re[i],xt.im[i]);
    //fprintf(fp2,"\n %d  %le   %le ", i,xt.re[i], xt.im[i]);
    fprintf(fp2,"\n %d  %le   ", i,xt.re[i] );
    }
    printf ("\n   input   '0' for exit ");
    std::cin>>key;
    return 0;

    }

    // stdafx.cpp : source file that includes just the standard includes
    // synthsig.pch will be the pre-compiled header
    // stdafx.obj will contain the pre-compiled type information

    #include "stdafx.h"

    // TODO: reference any additional headers you need in STDAFX.H
    // and not in this file

    // stdafx.h : include file for standard system include files,
    // or project specific include files that are used frequently, but
    // are changed infrequently
    //

    #pragma once


    #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
    #include <stdio.h>
    #include <tchar.h>

    // TODO: reference additional headers your program requires here

    spectrum.txt

      44100    16

     0   0   0
     1   0   0
     2   10   0
     3   0   0
     4   0   0
     5   0   0
     6   0   0
     7   0   0
     8   0   0
     9   0   0
     10   0   0
     11   0   0
     12   0   0
     13   0   0
     14  0   0
     15  0   0

    Friday, January 26, 2018 1:27 PM
  • // orthosynth.cpp : Defines the entry point for the console application.
    //

    #include "stdafx.h"


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

    #define NMAX 1000
    #define M_PI 3.1415926535897932384626433832795

    typedef struct {
    double re[NMAX];
    double im[NMAX];
    } CData;

     typedef struct  {
         double data[NMAX];
     } Array;


    CData  cdft(CData xn  ,int N )
     {
       
     
     CData Xk ;
     int k,n  ;
     double angle;

     for (k=0;k<N;k++)
      {
      Xk.re[k]=0;
      Xk.im[k]=0;
      for(n=0;n<N;n++)
        {
        angle=-2*M_PI*k*n/N;  
        Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
        Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
        }

     }

     return Xk;
     }

     CData  cidft(CData idft_Xk, int N  )
     {

     int k,n ;
     double angle;
     CData  idft_xn;
     
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (n=0;n<N;n++)
      {
      idft_xn.re[n]=0;
      idft_xn.im[n]=0;
      for(k=0;k<N;k++)
       {
     
       angle=2*M_PI*k*n/N;
       idft_xn.re[n]=idft_xn.re[n]+idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
       idft_xn.im[n]=idft_xn.im[n]+idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
       }

     idft_xn.im[n]=idft_xn.im[n]/N;
     idft_xn.re[n]=idft_xn.re[n]/N;

     }

    return idft_xn ;
     }

    CData ScaleDFT(CData xk,int  num)
     {

     int i;
     for (i=0; i<num ; i++)
       {
        xk.re[i]=xk.re[i]/num;
        xk.im[i]=xk.im[i]/num;
       }
     return xk;
     }

    CData ScaleIQidft(CData Xk, int num)
     {
     int i;
       for (i=0;  i<num ; i++)
      {
       Xk.re[i]=Xk.re[i]*num;
       Xk.im[i]=Xk.im[i]*num;
      }

    return Xk;

    }

     CData hilbert(double *xr , int N)
     {
      int  k,n,i, no2 ;
      double angle;
      double h[NMAX];
      double dft__xn_re[NMAX];
      double dft__xn_im[NMAX];
       double dft__Xk_re[NMAX];
       double dft__Xk_im[NMAX];
         double idft__XK_re[NMAX];
         double idft__XK_im[NMAX];
       CData idft__XN;

     
      if (N==0)
      {
      printf("\n Error");
      idft__XN.re[0]=0;
      idft__XN.im[0]=0;
      return idft__XN;
      }
      //xnreal=real(xr);
      no2=(int)floor ((double)N/2);
     
      for (i=0;i<N;i++)
      {
       dft__xn_re[i]=xr[i];
       dft__xn_im[i]=0;
      }

     //fft
     // dft__Xk=dft(dft__xn,n);

    for (k=0;k<N;k++)
      {
      dft__Xk_re[k]=0;
      dft__Xk_im[k]=0;
      for(n=0;n<N;n++)
        {
        angle=-2*M_PI*k*n/N;
        dft__Xk_re[k]+= dft__xn_re[n]*cos(angle)- dft__xn_im[n]*sin(angle);
        dft__Xk_im[k]+= dft__xn_re[n]*sin(angle)+ dft__xn_im[n]*cos(angle);
        }

     }

     //X={Xk.re,Xk.im};

     for (i=0;i<=N;i++)
      {
      h[i]=0;
      }


      if((N>0)&&((2*floor((double)N/2))==(N))) // n is even
      {

      h[0]=1;
       h[N/2]=1;
       for(k=1;k<N/2;k++)
       {
       h[k]=2;
       }
      }
      else //n is odd
      {
        if(N>0)
        {
         h[0]=1;
         for(k=1;k<(N+1)/2;k++)
         {
         h[k]=2;
         }
        }
      }

     for(i=0;i<N;i++)
      {
      idft__XK_re[i]= dft__Xk_re[i]*h[i];
      idft__XK_im[i]= dft__Xk_im[i]*h[i];
      }


     //idft__XN= idft(idft__XK, n);
     
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (n=0;n<N;n++)
      {
      idft__XN.re[n]=0;
      idft__XN.im[n]=0;
      for(k=0;k<N;k++)
       {
       //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
       angle=2*M_PI*k*n/N;
       idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
       idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
       }
     idft__XN.im[n]=idft__XN.im[n]/N;
     idft__XN.re[n]=idft__XN.re[n]/N;

     }


    return idft__XN;
     }

    CData getanalyticsignal(double *xn, int num)
     {
     CData x_n,Hilb;
     int i;
       //x=x+jy,y=0
       //hilb(x)=a+jb
       //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
       //I=x(t), Q=j*Hilb(xre(t))

    Hilb=hilbert(xn,num);
     
     for (i=0; i<num; i++)
     {
      x_n.re[i]=xn[i];
      x_n.im[i]=Hilb.im[i];
     }
     
     return x_n;
     }


     CData GetIQfromR_fi(CData Xinp, int numb)
     {
     int i;
     CData res;
      for (i=0;i<numb;i++)
      {
       Xinp.im[i]=Xinp.im[i]*M_PI/180;
      res.re[i]= Xinp.re[i]*cos(Xinp.im[i]);
      res.im[i]= Xinp.re[i]*sin(Xinp.im[i]);
      }

     return res;
     }

    double Arg(double re, double im)
    {
     
    return atan2(im,re); 
    }

    CData  get_mod_fi(CData x, int n, int mode   )
     {
      int i;
      CData res;
      double fi;
      for (i=0;i<n;i++)
      {
      res.re[i]=pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
      //double treshh=1e9;
      // x.im[i]=x.im[i]*treshh;
      // x.im[i]=floor(x.im[i])/treshh;
      // x.re[i]=x.re[i]*treshh;
      // x.re[i]=floor(x.re[i])/treshh;
      res.im[i]=Arg(x.re[i],x.im[i]);
      if (mode==1)  { res.im[i]=res.im[i]*180/M_PI;  }
              
      }
      return  res;
     }

    CData GetSpectrum(double *xn, int num)
    {
     
    CData xn_iq, Xf_IQ, Xf_modfi;
    xn_iq=getanalyticsignal( xn,  num);
     
    Xf_IQ=cdft(xn_iq  ,num );
    Xf_IQ=ScaleDFT(Xf_IQ,num);
    Xf_modfi=get_mod_fi(Xf_IQ, num,1   );
     
    return Xf_modfi;
    }

    CData SynthSignal(CData Xf_modfi, int num)

    {
    CData Xf_IQ,xn_iq;
    Xf_IQ=GetIQfromR_fi( Xf_modfi, num);
    xn_iq=cidft(Xf_IQ,num  );
    xn_iq=ScaleIQidft(xn_iq,  num);
    return xn_iq;
    }

    Array getenvelope(double *x, int num)
     {
     int i;
     CData  y;
     Array env;
     y=hilbert(x,num);
     for (i=0;i<num;i++)
      {
      env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
      }
      return env;
     }

    FILE   *fp1,*fp2 ;

    //int _tmain(int argc, _TCHAR* argv[])
     int main ()
    {

    double xin[NMAX];
    CData  Xf;
    CData xt;

    int k,n,i,N1;
    double fs=44100,T,x,f;
    int  key=0;

    printf ("\n For input  x(t) from file input '1' ");
    printf ("\n For manual  input  x(t) input   '0' ");
    std::cin>>key;
    switch(key)
    {
    case  1 :
         goto label_05;
      break;
    case  0 :
         goto label_06;
      break;
      
    }

    label_05:

     printf("\n Open spectrum.txt for reading Xf   ");
    if ((fp1=fopen("spectrum.txt","r"))==NULL)  printf ("\n File opening error ");
    fscanf(fp1,"%le",&fs);
    fscanf(fp1,"%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     fscanf (fp1,"%d",&k);
     fscanf (fp1,"%le", &Xf.re[n]);
     fscanf (fp1,"%le", &Xf.im[n]);
    }


    fclose(fp1);

    goto label_07;

    label_06:

    printf("\n Input fs,Hz , ( %lf ) ",fs) ;
    scanf( "%le",&fs);
    printf("\n Input number of samples of xin[n]  ") ;
    scanf( "%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     
     f=fs*n/N1; 
     printf ("\n n=%d ; f=%le Hz; ",n,f);
     printf ("\n Input |X(jw)| :" );
     scanf ( "%le",&Xf.re[n]);
     printf ("\nInput arg(X(jw)), deg :" );
     scanf ( "%le",&Xf.im[n]);
    }

     label_07:
    printf( "\n fs=%le Hz ",fs);
    printf( "\n Nsamples= %d ",N1);

    for (n=0 ; n<N1;n++)
    {
      f=fs*n/N1;
     printf ("\n f=%le Hz;|Xf(jw)|=%le;arg(Xf)=%le deg", f,Xf.re[n],Xf.im[n]);
    }

    printf("\n Input 0 for continue ");
    scanf("%s",&key);
    printf("\n OK.   ");


    /*****************************/

    label1:
     
    printf("\n Parsing   Xf->x[n]");
    xt=SynthSignal( Xf , N1);

    printf("\n Signal (xinput.txt): ");
     
     
    if ((fp2=fopen("orthoxn.txt","a"))==NULL)  printf ("\n File opening error");
    fprintf(fp2,"\n   %le   %d \n",fs,N1);

    for (i=0; i<N1;i++)
    {
     
    printf ("\n  %d ;   %le ;   %le  ",  i,xt.re[i],xt.im[i]);
     fprintf(fp2,"\n %d  %le   %le ", i,xt.re[i], xt.im[i]);
    //fprintf(fp2,"\n %d  %le   ", i,xt.re[i] );
    }
    printf ("\n   input   '0' for exit ");
    char key1;
    std::cin>>key1;
    std::cout << key1;


    return 0;

    }

    // stdafx.cpp : source file that includes just the standard includes
    // orthosynth.pch will be the pre-compiled header
    // stdafx.obj will contain the pre-compiled type information

    #include "stdafx.h"

    // TODO: reference any additional headers you need in STDAFX.H
    // and not in this file

    // stdafx.h : include file for standard system include files,
    // or project specific include files that are used frequently, but
    // are changed infrequently
    //

    #pragma once


    #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
    #include <stdio.h>
    #include <tchar.h>

    // TODO: reference additional headers your program requires here

    Friday, January 26, 2018 1:52 PM
  • // spectrumratio.cpp : Defines the entry point for the console application.
    //

    #include "stdafx.h"

     

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

    #define NMAX 1000
    #define M_PI 3.1415926535897932384626433832795

    typedef struct {
    double re[NMAX];
    double im[NMAX];
    } CData;

     typedef struct  {
         double data[NMAX];
     } Array;


    CData  cdft(CData xn  ,int N )
     {
       
     
     CData Xk ;
     int k,n  ;
     double angle;

     for (k=0;k<N;k++)
      {
      Xk.re[k]=0;
      Xk.im[k]=0;
      for(n=0;n<N;n++)
        {
        angle=-2*M_PI*k*n/N;  
        Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
        Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
        }

     }

     return Xk;
     }

     CData  cidft(CData idft_Xk, int N  )
     {

     int k,n ;
     double angle;
     CData  idft_xn;
     
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (n=0;n<N;n++)
      {
      idft_xn.re[n]=0;
      idft_xn.im[n]=0;
      for(k=0;k<N;k++)
       {
     
       angle=2*M_PI*k*n/N;
       idft_xn.re[n]=idft_xn.re[n]+idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
       idft_xn.im[n]=idft_xn.im[n]+idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
       }

     idft_xn.im[n]=idft_xn.im[n]/N;
     idft_xn.re[n]=idft_xn.re[n]/N;

     }

    return idft_xn ;
     }

    CData ScaleDFT(CData xk,int  num)
     {

     int i;
     for (i=0; i<num ; i++)
       {
        xk.re[i]=xk.re[i]/num;
        xk.im[i]=xk.im[i]/num;
       }
     return xk;
     }

    CData ScaleIQidft(CData Xk, int num)
     {
     int i;
       for (i=0;  i<num ; i++)
      {
       Xk.re[i]=Xk.re[i]*num;
       Xk.im[i]=Xk.im[i]*num;
      }

    return Xk;

    }

     CData hilbert(double *xr , int N)
     {
      int  k,n,i, no2 ;
      double angle;
      double h[NMAX];
      double dft__xn_re[NMAX];
      double dft__xn_im[NMAX];
       double dft__Xk_re[NMAX];
       double dft__Xk_im[NMAX];
         double idft__XK_re[NMAX];
         double idft__XK_im[NMAX];
       CData idft__XN;

     
      if (N==0)
      {
      printf("\n Error");
      idft__XN.re[0]=0;
      idft__XN.im[0]=0;
      return idft__XN;
      }
      //xnreal=real(xr);
      no2=(int)floor ((double)N/2);
     
      for (i=0;i<N;i++)
      {
       dft__xn_re[i]=xr[i];
       dft__xn_im[i]=0;
      }

     //fft
     // dft__Xk=dft(dft__xn,n);

    for (k=0;k<N;k++)
      {
      dft__Xk_re[k]=0;
      dft__Xk_im[k]=0;
      for(n=0;n<N;n++)
        {
        angle=-2*M_PI*k*n/N;
        dft__Xk_re[k]+= dft__xn_re[n]*cos(angle)- dft__xn_im[n]*sin(angle);
        dft__Xk_im[k]+= dft__xn_re[n]*sin(angle)+ dft__xn_im[n]*cos(angle);
        }

     }

     //X={Xk.re,Xk.im};

     for (i=0;i<=N;i++)
      {
      h[i]=0;
      }


      if((N>0)&&((2*floor((double)N/2))==(N))) // n is even
      {

      h[0]=1;
       h[N/2]=1;
       for(k=1;k<N/2;k++)
       {
       h[k]=2;
       }
      }
      else //n is odd
      {
        if(N>0)
        {
         h[0]=1;
         for(k=1;k<(N+1)/2;k++)
         {
         h[k]=2;
         }
        }
      }

     for(i=0;i<N;i++)
      {
      idft__XK_re[i]= dft__Xk_re[i]*h[i];
      idft__XK_im[i]= dft__Xk_im[i]*h[i];
      }


     //idft__XN= idft(idft__XK, n);
     
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (n=0;n<N;n++)
      {
      idft__XN.re[n]=0;
      idft__XN.im[n]=0;
      for(k=0;k<N;k++)
       {
       //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
       angle=2*M_PI*k*n/N;
       idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
       idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
       }
     idft__XN.im[n]=idft__XN.im[n]/N;
     idft__XN.re[n]=idft__XN.re[n]/N;

     }


    return idft__XN;
     }

    CData getanalyticsignal(double *xn, int num)
     {
     CData x_n,Hilb;
     int i;
       //x=x+jy,y=0
       //hilb(x)=a+jb
       //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
       //I=x(t), Q=j*Hilb(xre(t))

    Hilb=hilbert(xn,num);
     
     for (i=0; i<num; i++)
     {
      x_n.re[i]=xn[i];
      x_n.im[i]=Hilb.im[i];
     }
     
     return x_n;
     }


     
    double Arg(double re, double im)
    {
     
    return atan2 (im,re)  ; 
    }

    CData  get_mod_fi(CData x, int n, int mode   )
     {
      int i;
      CData res;
      double fi;
      for (i=0;i<n;i++)
      {
      res.re[i]=pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
      //double treshh=1e9;
      // x.im[i]=x.im[i]*treshh;
      // x.im[i]=floor(x.im[i])/treshh;
      // x.re[i]=x.re[i]*treshh;
      // x.re[i]=floor(x.re[i])/treshh;
      res.im[i]=Arg(x.re[i],x.im[i]);
      if (mode==1)  { res.im[i]=res.im[i]*180/M_PI;  }
              
      }
      return  res;
     }

    CData GetSpectrum(double *xn, int num)
    {
     
    CData xn_iq, Xf_IQ, Xf_modfi;
    xn_iq=getanalyticsignal( xn,  num);
     
    Xf_IQ=cdft(xn_iq  ,num );
    Xf_IQ=ScaleDFT(Xf_IQ,num);
    Xf_modfi=get_mod_fi(Xf_IQ, num,1   );
     
    return Xf_modfi;
    }

    Array getenvelope(double *x, int num)
     {
     int i;
     CData  y;
     Array env;
     y=hilbert(x,num);
     for (i=0;i<num;i++)
      {
      env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
      }
      return env;
     }
     
     
     
    double getmax( Array   k , int N)
    {
    double kmax=k.data[0];
    double k1=0;
    int i;
    for (i=1;i<N;i++)
    {
     k1=k.data[i];
      if(fabs( k1 )>fabs(kmax)) { kmax= k.data[i]; } 
      
    }
     
     return kmax;
    }
     
    double getmin( Array  k, int N)
    {
    double kmin=k.data[0];
    int i;
    for (i=0;i<N;i++)
    {
      if(fabs(k.data[i])<fabs(kmin)) { kmin= k.data[i]; } 
    }
     return kmin;
    }

    double getmid( Array  k,int N)
    {
    double kmid=0;
    int i;
    for (i=0;i<N;i++)
    {
     kmid+=   k.data[i];
    }
    kmid=kmid/N;
     return kmid;
    }


    FILE   *fp1,*fp2 ;

    int _tmain(int argc, _TCHAR* argv[])
    //int main ()
    {
    double delta =1e-15;
    double xin[NMAX];
    double yout[NMAX];
    CData Xf1,Xf2,Xf3;
    double Kmax=0, Kmin=0, Kmid=0;
    int key;
    int k,n,i,N1,N2;
    double fs1,fs2,T,x,f;


    /***********input x[n] *************/
    label_00:
    printf ("\n For input  x[n] from file input '1' ");
    printf ("\n For manual  input  x[n] input   '0' ");
     std::cin>>key;
    if (key== 1  ) { goto label_01;}
    if (key== 0  ) { goto label_02;}

    label_01:

     printf("\n Open xinput.txt for reading x[n]   ");
    if ((fp1=fopen("xinput.txt","r"))==NULL)  printf ("\n File opening error ");
    fscanf(fp1,"%le",&fs1);
    fscanf(fp1,"%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     fscanf (fp1,"%d",&k);
     fscanf (fp1,"%le",&xin[n]);
    }


    fclose(fp1);

    goto label_03;

    label_02:

    printf("\n Input fs input,Hz  " ) ;
    scanf( "%le",&fs1);
    printf("\n Input number of samples of xin[n]  ") ;
    scanf( "%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     T=n*1/(N1*fs1); 
     printf ("\n n=%d ; t=%le s; xin[n]= ",n,T);
     scanf ( "%le",&xin[n]);
    }

     label_03:
    printf( "\n fs1=%le Hz ",fs1);
    printf( "\n Nsamples1= %d ",N1);

    for (n=0 ; n<N1;n++)
    {
     T=n*1/(N1*fs1); 
     printf("\n n=%d x[n]=%le,  t=%le s",n,xin[n],T);
    }

    printf("\n Input 0 for continue ");
    char key1;
    std::cin>>key1;
    printf("\n OK.   ");

    /***********input y[n] *************/

    printf ("\n For input  y[n] from file input '1' ");
    printf ("\n For manual  input  y[n] input   '0' ");
    std::cin>>key;
    if (key== 1  ) { goto label_04;}
    if (key== 0  ) { goto label_05;}

    label_04:

     printf("\n Open youtput.txt for reading y[n]   ");
    if ((fp1=fopen("youtput.txt","r"))==NULL)  printf ("\n File opening error ");
    fscanf(fp1,"%le",&fs2);
    fscanf(fp1,"%d", &N2);

    for (n=0 ; n<N2;n++)
    {
     fscanf (fp1,"%d",&k);
     fscanf (fp1,"%le",&yout[n]);
    }


    fclose(fp1);

    goto label_06;

    label_05:

    printf("\n Input fs out,Hz , ( %lf ) ",fs1) ;
    scanf( "%le",&fs2);
    printf("\n Input number of samples of yout[n]  ") ;
    scanf( "%d",&N2);

    for (n=0 ; n<N2;n++)
    {
     T=n*1/(N2*fs2); 
     printf ("\n n=%d ; t=%le s; yout[n]= ",n,T);
     scanf ( "%le",&yout[n]);
    }

     label_06:


    if (N1!=N2) printf("Nin <> Nout");

    printf( "\n fs =%le Hz ",fs2);
    printf( "\n Nsamples  out = %d ",N2);

    for (n=0 ; n<N2;n++)
    {
     T=n*1/(N2*fs2); 
     printf("\n n=%d y[n]=%le,  t=%le s",n,yout[n],T);
    }

    printf("\n Input 0 for continue ");
    std::cin>>key;
    printf("\n OK.   ");

    /***********************/
    if ((fs2!=fs1)||(N1!=N2)) { printf("Nin <> Nout"); goto label_00;}


    /*****************************/

    label1:
     
    printf("\n Obtaining spectrum of x[n]");
    Xf1=GetSpectrum(xin, N1) ;
    printf("\n Obtaining spectrum of y[n]");
    Xf2=GetSpectrum(yout, N2) ;


    printf("\n  |Sy(jw)|/|Sx(jw)| : ");
     
    if ((fp2=fopen("ratio.txt","a"))==NULL)  printf ("\n File opening error");
    fprintf(fp2,"\n\n   %le   %d \n",fs1,N1);


    for (i=0; i<N1;i++)
    {
      if (Xf1.re[i]==0 ) { Xf3.re[i]=Xf2.re[i]/(Xf1.re[i]+delta); }
      if (Xf1.re[i]!=0 ) { Xf3.re[i]=Xf2.re[i]/(Xf1.re[i] ); }
     Xf3.im[i]=Xf2.im[i]-Xf1.im[i];

     
    printf ("\n %d ;|Sy|/|Sx|=%le ;arg(Sy(jw))-arg(Sx(jw))=%le deg ",  i,Xf3.re[i],Xf3.im[i]);
    fprintf (fp2,"\n %d   %le  %le",  i,Xf3.re[i], Xf3.im[i]);

    }
    fclose(fp2);

     

    printf ("\n  input   '0' for exit ");
    char key2;
    std::cin>>key2;

    return 0;

    }

    // stdafx.cpp : source file that includes just the standard includes
    // spectrumratio.pch will be the pre-compiled header
    // stdafx.obj will contain the pre-compiled type information

    #include "stdafx.h"

    // TODO: reference any additional headers you need in STDAFX.H
    // and not in this file

    // stdafx.h : include file for standard system include files,
    // or project specific include files that are used frequently, but
    // are changed infrequently
    //

    #pragma once


    #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
    #include <stdio.h>
    #include <tchar.h>

    // TODO: reference additional headers your program requires here

    Friday, January 26, 2018 2:04 PM
  • make.bat

    set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"

     
     javac THDObtainer.java
     jar cvfe THDObtainer.jar THDObtainer   *.class
     jar uf  THDObtainer.jar THDObtainer.java
    pause 0

     

    set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"
     rem java  THDObtainer
     java -jar THDObtainer.jar

    pause 0

    import java.util.Scanner;
    import java.io.*;


       class  THDParser
     {

     

     
     
    static int  NMAX = 1024;

      public static double[] Getfassoc(double fs , int  Numb )
    {
       double[] fvect  =new double[Numb];
      
       for(int k = 0;k<Numb;k++)
      {
         fvect[k] = fs * k / Numb;
     
      }
     return fvect;
    }
     
     
    public static double[] Gettassoc(double  fs ,    int Numb )
     {
      double[] tvect=new double[Numb];
       for(int N = 0 ; N<Numb; N++)
        {
         tvect[N] = N/fs ;
        }
     return tvect;
    }


      public static double  Arg(double x_re  , double x_im  )
      {
      double Arg=0;
     if (x_re > 0)  { Arg = Math.atan(x_im / x_re); }
     if ((x_re == 0) &&  (x_im > 0)) { Arg = Math.PI / 2;}
     if ((x_re == 0) &&  (x_im < 0)) { Arg = Math.PI / 2 ; }
     if ((x_re < 0) &&  (x_im >= 0)) { Arg = Math.atan(x_im / x_re) + Math.PI; }
     if ((x_re < 0) &&  (x_im < 0))  { Arg = Math.atan(x_im / x_re) - Math.PI; }
     if ((x_re == 0) &&  (x_im == 0))  {
                                      Arg = 0;
                                   /* MsgBox (" Xre=0 and Xim=0, Invalid result"); */
                                    }
     
    return Arg;
    }


     /*
     '36.  - 4. + 9.6568542i  - 4. + 4.i  - 4. + 1.6568542i  - 4.  - 4. - 1.6568542i
     
      '         column 7 to 8
      
      ' - 4. - 4.i  - 4. - 9.6568542i
    */

      public static double[][]  GetIQfromR_fi(  int Num  ,double [][] Xk, int mode )
    {
       //long int i ;
        
          double arg1 ;
          double[][] Xk1=new double [Num][2];
        arg1 = 1;
        if( mode == 1) { arg1 = Math.PI / 180;}       
        for(int i = 0 ;i< Num;i++)
         {
          Xk1[i][0] =  Xk[i][0]*Math.cos(arg1*Xk[i][1]);
          Xk1[i][1] =  Xk[i][0]*Math.sin(arg1*Xk[i][1]);
         }
     return Xk;   
    }

     


      public static double[][] CDFT(int numdft ,double [][] dft_xn )
    {
      int i, k, N;
     
     double angle ;
     double[][] dft_Xk=new double [numdft ][2];
    for( k = 0 ;k<numdft ;k++)
     {
      dft_Xk[k][0] = 0;
      dft_Xk[k][1] = 0;
       for( N = 0;N<numdft;N++ )
       {
         angle = -2*Math.PI*k*N/numdft;
       dft_Xk[k][0]+=dft_xn[N][0]*Math.cos(angle)-dft_xn[N][1]*Math.sin(angle);
       dft_Xk[k][1]+=dft_xn[N][0]*Math.sin(angle)+dft_xn[N][1]*Math.cos(angle);
       }
      }
    return dft_Xk;
    }
     
    public static double[][]  ScaleDFT(int Num , double [][]  Xk )
    {
     
      for(int i = 0 ;i<Num;i++)
      {
         Xk[i][0] = Xk[i][0]/Num;
         Xk[i][1] = Xk[i][1]/Num;
      }
     return Xk;
    }

    public static double[][] ScaleIQidft(int Num  , double[][] ScXk)
    {
      
       for (int i = 0;i<Num;i++)
       {
        ScXk[i][0] = ScXk[i][0]*Num;
        ScXk[i][1] = ScXk[i][1]*Num;
       }

    return ScXk;
    }

    public static double[][] CIDFT(int numidft , double[][] idft_Xk)
    {
      double angle;
      double[][] idft_xn=new double[numidft][2] ;
      for(int n = 0 ; n<numidft; n++)
      {
     
       idft_xn[n][0] = 0;
       idft_xn[n][1] = 0;
       for (int k=0; k<numidft; k++)
       {
        angle =2*Math.PI*k*n/numidft;
        idft_xn[n][0]+=idft_Xk[k][0]*Math.cos(angle)-idft_Xk[k][1]*Math.sin(angle);
        idft_xn[n][1]+=idft_Xk[k][0]*Math.sin(angle)+idft_Xk[k][1]*Math.cos(angle);
       }
       
        idft_xn[n][0]= idft_xn[n][0]/numidft;
        idft_xn[n][1]= idft_xn[n][1]/numidft;
     }
     
     return idft_xn;
    }

    public  static double[][] Hilbert(double[] xre , int numb  )
    {
      /* MATLAB algorithm */
       int k, i, no2;
       double[] h= new double [NMAX+1];
     
     
       double[][] dft_xn   =new  double [NMAX][2];
       double[][] dft_Xk   =new  double [NMAX][2];
       double[][] idft_Xk  =new  double [NMAX][2];
       double[][] idft_xn  =new  double [NMAX][2];

      if (numb == 0)
       {
               idft_xn[0][0]=0;
               idft_xn[0][1]=0;
        
       /*MsgBox (" Hilbert() :   Parsing error ")*/
        }


       no2 = (int) Math.floor(numb / 2); //
      
       for(i = 0;i<numb; i++)
       {
        dft_xn[i][0]=xre[i]; 
        dft_xn[i][1]=0; 
       }


       dft_Xk=CDFT(numb, dft_xn ); 
       i=0;
       for (i=0;i <= numb ;i++){ h[i]=0; }

      /*if  N is even*/
     
        if( (numb> 0)&&((2*no2)==numb) )
        {
            h[0] = 1; // 0
            h[no2] = 1;// f Niquist
             for( k = 1; k<no2;k++)
             {
                h[k]= 2;
             }
        }

       /* if N is odd */

        else {
             if(numb > 0)
              {

                h[0]  = 1;
                for( k = 1; k<=((0.5*(numb + 1))-1); k++)
                {
                 h[k] = 2;
                }
           
               }
              }


         for( i=0;i<numb;i++)
          {
           idft_Xk[i][0]=dft_Xk[i][0]*h[i];
           idft_Xk[i][1]=dft_Xk[i][1]*h[i];
          }
     
        
         idft_xn=CIDFT(numb, idft_Xk);
       
     return idft_xn;
     
    }


    public static double[][] GetAnalyticSignal(double [] xnre , int Numb )
    {
     //  'I=x(t)
       // 'Q=jHilb(t)
     //  'Hilb=Hilbert()
       
       double[][] Hilb= new double[Numb][2];
       double[][] x_iq= new double[Numb][2];

       Hilb= Hilbert(xnre, Numb );
       for(int i = 0; i<Numb;i++)
       {
        x_iq[i][0] = xnre[i] ;
        x_iq[i][1] = Hilb[i][1];
       }
    return  x_iq;
    }

    public static double[] getenvelope(double [] xnre , int Numb)
       {
      
       
       double[][] Hilb= new double[Numb][2];
       double[]  env= new double[Numb] ;
       Hilb= Hilbert(xnre, Numb );

       for(int i = 0 ; i<Numb; i++)
       {
       env[i]=Math.sqrt((Hilb[i][0]*Hilb[i][0])+(Hilb[i][1]*Hilb[i][1]));
       }
    return env;

    }


    public static double[][] GetSpectrumIQ(double[] xn , int Num  )
     {
      double[][] x_iq= new double[Num ][2];
      double[][] Xf_iq= new double[Num ][2];
      x_iq=GetAnalyticSignal(xn , Num  ); 
      Xf_iq=CDFT(Num, x_iq);
      Xf_iq=ScaleDFT(Num,Xf_iq);
    return Xf_iq;
    }


    public static double[][] GetR_fi(int  Numb  , double[][] Xf_iq , int mode )
    {

        double[][] Xf_rfi=new double [Numb][2];
        for (int i = 0 ;i<Numb;i++)
       {
        Xf_rfi[i][0] =Math.sqrt((Xf_iq[i][0]*Xf_iq[i][0])+(Xf_iq[i][1]*Xf_iq[i][1]) );
        Xf_rfi[i][1] = Arg(Xf_iq[i][0], Xf_iq[i][1]);
        if( mode == 1) { Xf_rfi[i][1] = Xf_rfi[i][1]*180/Math.PI; }
       }

    return Xf_rfi;
    }


    public static double[][] GetSpectrum(double[] xn , int Num ,  int mode  )
    {
              double[][] x_iq= new double[Num][2];
              double[][] Xf_iq= new double[Num][2];
              double[][] Xf_rfi=new double [Num][2];

           x_iq=GetAnalyticSignal(xn, Num );
           Xf_iq=CDFT(Num,x_iq);
           Xf_iq=ScaleDFT(Num, Xf_iq);
           Xf_rfi=GetR_fi(Num, Xf_iq, mode);
      return Xf_rfi;
    }
     
     
    public static double[][] GetSpectrumfromIQ(double[][] x_iq,int Num, int mode)
    {  
           double[][] Xf_iq= new double[Num][2];
           double[][] Xf_rfi=new double [Num][2];

          Xf_iq=CDFT(Num,x_iq);
          Xf_iq=ScaleDFT(Num, Xf_iq);
          Xf_rfi=GetR_fi(Num,Xf_rfi, mode);
      return Xf_rfi;
    }


    public static double[][] SynthSignalfromIQ(int Num  , double[][] Xk_iq)
    {
        double[][] xn_iq= new double[Num][2];
        
              xn_iq= CIDFT(Num, Xk_iq);
              xn_iq= ScaleIQidft(Num, xn_iq);
    return xn_iq;
    }


    public static double[][] SynthSignalfromRFi( double[][] Xk_rfi, int Num, int mode)
    {
    double[][] Xk_iq= new double [Num][2];
    Xk_iq=GetIQfromR_fi(Num, Xk_rfi, mode );
    return SynthSignalfromIQ(Num, Xk_iq);
    }

    public static void printcomplex(double[][] xn  , int num, int col1,int col2)

        
         for (int i=0 ;i<num;i++)
            {
             
              System.out.print("\n i= ");
              System.out.print(i);
              System.out.print(" ");
              System.out.print(String.format("%.08f",xn[i][0]));
              System.out.print(" ");
              System.out.print(String.format("%.08f",xn[i][1]));
        
            }
     System.out.print("\n ");
    return;
    }

    public static void printarray(double[] xn , int  Num  , int col1   )
    {
         for (int i = 0; i<Num;i++)
          {
            System.out.print("\n ");
            System.out.print(String.format("%.08f",xn[i]));
          
          }
     System.out.print("\n ");
    return;
    }

    public static double deltat(double[] argF , double f0 , double fs , int N )
     {
      double[] fv =new double[N+1];
      int  N0;
      double deltat;
       fv=Getfassoc(fs, N);
       N0 = 0;
     for(int i = 0;i<N-1;i++)
     {
     if ((fv[i]>= f0)&&(fv[i+1]<=f0)) { N0= i; }
     }
     deltat =-(argF[N0] -argF[0])/(2*Math.PI*fv[N0]-2*Math.PI*fv[0]);
    return deltat;
    }


     
    public static void writetofile(double[][] xnc  ,  String FilePath , int N ) throws IOException
    {  
     File outFile1 = new File (FilePath);
     FileWriter fWriter1 = new FileWriter (outFile1);
     PrintWriter pWriter1 = new PrintWriter (fWriter1);

        for(int i = 0 ;i< N;i++)
        {
             pWriter1.println("\n ");
             pWriter1.println("  ");
      pWriter1.print(String.format("%d",i) );
             pWriter1.print("   "); 
      pWriter1.print(String.format("%e",xnc[i][0] ));
      pWriter1.print("   ");
             pWriter1.print(String.format("%e",xnc[i][1] ));
      pWriter1.print("   ");
       }
      pWriter1.close();
     return;        
      
     }

     


    public static void TestHilbert( )
    {

     int i ;
     int N = 8;
     double  fs  ;
     double[] fv=new double[N ];
     double[] tv=new double[N ];
     double[][] Xf=new double[N ][2];

     double[][] Yf=new double[N ][2];
     double[][] Hilb=new double[N][2];
     double[][] Idft=new double[N ][2];
     double[][] SpX=new double[N][2];
     
     double[] env=new double[N];
     double[][] Xn=new double[N][2];
     double[] Xn1=new double[N];
     
     double[] x_re=new double[N];
     double[][] Yn=new double[N][2];
     
     
     Xf[0][0] = 1;
     Xf[1][0] = 2;
     Xf[2][0] = 3;
     Xf[3][0] = 4;
     Xf[4][0] = 5;
     Xf[5][0] = 6;
     Xf[6][0] = 7;
     Xf[7][0] = 8;

     Xf[0][1] = 0;
     Xf[1][1] = 0;
     Xf[2][1] = 0;
     Xf[3][1] = 0;
     Xf[4][1] = 0;
     Xf[5][1] = 0;
     Xf[6][1] = 0;
     Xf[7][1] = 0;
     
      fs = 1000;

      printcomplex(Xf,N,1,2);
       fv=Getfassoc(fs, N);
       tv=Gettassoc(fs, N);
    System.out.print("\n f= ");
     printarray(fv, N,  0);
    System.out.print("\n t= ");
     printarray(tv, N,  0);


      Yn=CDFT(N,Xf);
    System.out.print("\n CDFT(XF)= ");
      printcomplex(Yn, N, 5, 6);
      Idft=CIDFT(N, Yn );
    System.out.print("\n CIDFT(Yn)= ");
      printcomplex(Idft, N, 7, 8);
        for(i=0;i<N;i++)
      {
      Xn1[i]=Idft[i][0];
      }

       Hilb=Hilbert(Xn1, N );
    System.out.print("\n Hilbert(Xn1)= ");
       printcomplex(Hilb, N, 9, 10);
       SpX=GetSpectrum(Xn1, N , 1);
    System.out.print("\n Get spectrum= ");
       printcomplex(SpX, N, 11, 12);
      
     System.out.print("\n Get envelope (Xn1) ");
       env=getenvelope(Xn1, N);
       printarray(env, N, 13);
      
       System.out.print("\nPoke new  aray ");

      //  cos (90° – x) = sin x
      // -sin ( x-90) = cos x 
      // Hilbert(sin(x)) = -cos(t)
      // Hilbert(Cos(x)) = sin(x)
     
     // For i = 0 To n - 1
     // Xi1(i + 1) = 1 * Sin(2 * M_PI * 15.625 * i / fs + 0) + 0.5 * Cos(2 * M_PI * 93.75 * i / fs + 0)
     // Next i
      
       N = 16;

      


      double [][] Xn2=new double[N][2] ;
      double [][] Xf2=new double[N][2] ;
       for (i= 0;i<N; i++)
        {
        Xf2[i][0] = 0;
        Xf2[i][1] = 0;
        }
       
       Xf2[1][0] = 1;
       Xf2[1][1] = -45;
       Xf2[5][0] = 0;
       Xf2[5][1] = 0;
       Xn2 =SynthSignalfromRFi(Xf2, N, 1);
      System.out.print("\n Get Xn from Xf");
      printcomplex(Xn2, N, 14, 15);
      fv=Getfassoc(fs, N);
     System.out.print("\n f[]=");
       printarray(fv, N, 16);
        System.out.print("\n SpX[]=");  
         SpX=GetSpectrumfromIQ(Xn2, N, 1);
         printcomplex(SpX, N, 17, 18);

     
     double[] Xn4=new double[N]; 
          for (i = 0 ; i<N; i++)
         {
          Xn4[i] = 0.5*Math.cos(2*Math.PI*250*i/fs+0);
          Xn4[i] +=1.5*Math.cos(2*Math.PI*125*i/fs+Math.PI/2);
          Xn4[i] +=2*Math.cos(2*Math.PI*375 *i/fs-Math.PI/2);
         }

      double [][] Xn3=new double[N][2] ;
      double [][] Xf3=new double[N][2] ;
      double [][] Yn3=new double[N][2] ;
     System.out.print("\n Spectrum (Xn4)=");    
         Xf3=GetSpectrum(Xn4, N , 1);
         printcomplex(Xf3, N, 19, 20);
      System.out.print("\n Spectrum IQ(Xn4)=");
         Yn3=GetSpectrumIQ(Xn4, N );
         printcomplex(Yn3, N, 21, 22);
     return;

    }


    public int  GetNumF1(double f1 ,int Numb  ,double fs )
    {
     return (int) (Math.floor(f1*Numb/fs));
    }

    public static double GetTHD_F(double[] xn , int Numb, double fs, double f1)
     {

    int numbf1;
    double[][] Xf =new double[Numb][2];
    double Urms1h;
    double Uhharm;
    int i ;

    //
    //f 0  5  10  15  20   25   30   35 Hz  N=8
    //i 0  1   2   3   4   5    6     7
    //         f1      2f1    3f1
    //      numb1=2
    //
    //f 0  5  10  15  20   25   30   35    40   45  50  55   60 Hz   N=13 imax=12
    //i 0  1   2   3   4   5    6     7    8    9   10  11   12
    //            f1           2f1             3f1           4f1
    //      numb1=3
    //

    numbf1=(int)Math.floor(f1*Numb/fs);//   2  , get associated vector of frequency
    Xf=GetSpectrum(xn, Numb,  0);
    Urms1h= Xf[numbf1][0]/Math.sqrt(2);
    i = numbf1*2;
     
    Uhharm=0;
    while (i<Numb)
    {
     if (i<Numb) { Uhharm+=(Xf[i][0]*Xf[i][0])/2 ; } 
     i=i+numbf1;
     
    }

    Uhharm =Math.sqrt(Uhharm) ;
     
    return Uhharm*100/(Urms1h);
    }


    public static double GetTHD_R(double[] xn , int Numb, double fs, double f1)
    {
    int numbf1;
    double[][]  Xf =new double[Numb][2];
    double Urmssigah;
    double Uhharm;
    int i ;

     
    numbf1=(int)Math.floor(f1*Numb/fs);  //   get associated vector of frequency, ORIGIN=0
    Xf=GetSpectrum(xn, Numb, 0);
    /* parser of Urms of all harmonics */

    // THD-R=Urmshharm*100/urmsig =sqrt( sum(Xrms[i] ^2))/Vrmssig;
    Urmssigah=0;
    i=numbf1;

     while (i<Numb)
     {
     if (i<Numb) { Urmssigah+=(Xf[i][0]*Xf[i][0])/2;  }
     i=i+1;
     }

    Urmssigah=Math.sqrt (Urmssigah );
    Uhharm=0;
     /* parser of Urms of higher harmonics */

     i=numbf1*2; // ' number of  bin of 2-nd harmonic 
      while (i<Numb)
     {
     if (i<Numb) { Uhharm+=(Xf[i][0]*Xf[i][0])/2;  }
     i=i+numbf1;         
     }

     Uhharm=Math.sqrt (Uhharm) ;
     
     return (Uhharm*100/Urmssigah);
    }

    public static double GetRMS(double[] xn , int Numb )
    {
    int i ;
    double Urms;
    Urms = 0;
    for (i=0;i<Numb;i++)
    {
     Urms+=xn[i]*xn[i];
    }
    Urms=Math.sqrt(Urms/Numb);

    return Urms;
    }


    public static double GetTHD_N(double[] xn , int Numb, double fs  , double f1 )
    {
    int numbf1;
    double[][] Xf=new double [Numb][2];
     
    double Urmssig;
    double Uhharmn;
    int i ;
    /*
    Urmssig = 0;
    for(i=0 ;i<Numb;i++)
    {
     Urmssig+=xn[i]*xn[i];
    }
    Urmssig=Math.sqrt(Urmssig/Numb);

    */
     numbf1=(int)Math.floor (f1*Numb/fs);
     Xf=GetSpectrum(xn, Numb, 0);

    i=0;
    Urmssig=0;
     while (i<Numb)
    {
    { Urmssig+=Xf[i][0]*Xf[i][0]/2;   }
    i=i+1;
    }

    Urmssig= Math.sqrt(Urmssig);

    Uhharmn=0;
    i=0;

     while (i<Numb)
    {
    if (i!= numbf1) { Uhharmn+=Xf[i][0]*Xf[i][0]/2;   }
    i=i+1;
    }

    Uhharmn= Math.sqrt(Uhharmn);
    return (Uhharmn*100/Urmssig);
    }

    public static void PrintValueDouble(String str, double s,String str1)
    {
     
    System.out.print(str);
    System.out.print(String.format("%.08f",s));
     System.out.print(str1);
    return;
    }

    public static void PrintValueInt(String str, int s,String str1)
    {
     
    System.out.print(str);
    System.out.print(String.format("%d",s));
     System.out.print(str1);
    return;
    }

    public static  double GetValueDouble(String s,Scanner scan1)
    {
     
    System.out.print(s);

    double result;
    result=(double)  scan1.nextDouble();
     
    return result;
    }


    public static int GetValueInt(String s,Scanner scan1)
    {
    System.out.print(s);
     
    int result=(int) scan1.nextInt();
     
    return result;
    }

    public static void TestTHDN( )
     {


     int N;
     Scanner sc2 = new Scanner (System.in);
    PrintValueInt("\n N=",16," ");
     N =16;// GetValueInt("\n Input N ",sc2);

     double [][] Xf =new double[N][2];
     double [][] X_iq=new  double[N ][2];
     double [][]  Xn =new  double[N ][2] ;
     double []  Xn1=new  double[N ] ;
     double [] fvect  =new double[N ] ;
     double fs;
     double f1;
     int i;

     fs =44100;// GetValueDouble("\n Input fs,Hz ",sc2);
     f1 =11025;//GetValueDouble("\n Input f [1-st harm] ,Hz ",sc2);
     fvect= Getfassoc(fs,  N);

     for(i = 0;i<N;i++)
     {
     
     //PrintValueDouble("\n f=",fvect[i]);
     Xf[i][0] =0;//GetValueDouble("\n input Um[f] " ,sc2);
     Xf[i][1] =0;//GetValueDouble("\n input fi[f] ",sc2);

     }
    PrintValueDouble("\n fs=",fs," ");
    PrintValueDouble("\n f1=",f1," ");

    PrintValueDouble("\n f=",fvect[4]," Hz");
    Xf[4][0]=1;
     PrintValueDouble("  |X[f]|=",Xf[4][0]," ");

    PrintValueDouble("\n f=",fvect[8]," Hz");
    Xf[8][0]=0.3;
     PrintValueDouble("  |X[f]|=",Xf[8][0]," ");

    PrintValueDouble("\n f=",fvect[9]," Hz");
    Xf[9][0]=0.1;
     PrintValueDouble("  |X[f]|=",Xf[9][0]," ");

     for(i = 0;i<N;i++)
     {
     
     PrintValueDouble("\n f=",fvect[i]," ");
     PrintValueDouble("  |X[f]|=",Xf[i][0], " ");
     PrintValueDouble("  fi=",Xf[i][1]," "); 

     }


      Xn=SynthSignalfromRFi(Xf,  N, 1);
      printcomplex(Xn, N, 5, 6);
     for (i=0;i<N;i++)
    {
    Xn1[i]=Xn[i][0];
    }
     
     double THD_N= GetTHD_N(Xn1, N, fs, f1);
     PrintValueDouble("\n THD+N= ", THD_N," %");
     double THD_R= GetTHD_R(Xn1, N, fs, f1);
      PrintValueDouble("\n THD-R= ",THD_R," %");
    double THD_F= GetTHD_F(Xn1, N, fs, f1 );
     PrintValueDouble("\n THD-F= ",THD_F," %");
    double RMS=GetRMS(Xn1, N);
     PrintValueDouble("\n U RMS=",RMS," units of ADC");
      sc2.close();
     return ;
     }

    }

    public   class  THDObtainer
     {
    public static void main(String[] args) throws IOException
     {
            THDParser.TestHilbert();
            THDParser.TestTHDN();
            }

    }

    Friday, January 26, 2018 2:08 PM
  • spectrumstat - utility for obtaining max. levels of signals on frequencies of the spectrum using short-time Fourier    transform and founder of the max. vaues

    xinputns.txt

    4,410000e+04   16  8
         
    0   1,000000e+01  
     
    1   -7,071068e+00  
     
    2   3,061617e-15  
     
    3   7,071068e+00  
     
    4   -1,000000e+01  
     
    5   7,071068e+00  
     
    6   -2,694842e-14  
     
    7   -7,071068e+00
    8   1,000000e+01  
     
    9   -7,071068e+00  
     
    10   3,061617e-15  
     
    11   7,071068e+00  
     
    12   -1,000000e+01  
     
    13   7,071068e+00  
     
    14   -2,694842e-14  
     
    15   -7,071068e+00   

    4,410000e+04   16

    maxspectr.txt
     
    0   3,020133e-15  0,000000e+00 
     
    1   1,330312e-07  5,512500e+03 
     
    2   6,447902e-15  1,102500e+04 
     
    3   1,000000e+01  1,653750e+04 
     
    4   4,642237e-15  2,205000e+04 
     
    5   5,151483e-15  2,756250e+04 
     
    6   2,334110e-15  3,307500e+04 
     
    7   7,772354e-15  3,858750e+04 

    make.bat

    set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"

     
     javac Spectrumstat.java
     jar cvfe Spectrumstat.jar Spectrumstat   *.class
      jar uf  Spectrumstat.jar Spectrumstat.java xinputns.txt
    pause 0

    run.bat

    set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"
     
     java -jar  Spectrumstat.jar  
    pause 0

    //javac Spectrum.java
    //java Spectrum

    import java.util.Scanner;
    import java.io.*;


     public class Spectrumstat
    {

          

      public static double[][]   cdft(double[][] xn  ,int N )
        {
      
       int j;
       double[][] Xk=new  double[N][2];
       int k,n  ;
       double angle;


     for (k=0;k<N;k++)
      {
      Xk[k][0]=0;
      Xk[k][1]=0;
      for(n=0;n<N;n++)
        {
        angle=-2*Math.PI*k*n/N;  
        Xk[k][0]+= xn[n][0]*Math.cos(angle)-xn[n][1]*Math.sin(angle);
        Xk[k][1]+= xn[n][0]*Math.sin(angle)+xn[n][1]*Math.cos(angle);
        }
     }

     return Xk;
     }


     public static double[][] ScaleDFT(double[][] xk,int  num)
      {
     
      for (int i=0; i<num ; i++)
       {
        xk[i][0] =xk[i][0] /num;
        xk[i][1] =xk[i][1] /num;
       }
     return xk;
     }

     public static double[][] hilbert(double[] xr  )
          {
          
          double angle;
          double[]  h=new double[xr.length+1]  ;
          double[]  dftxn_re=new double[xr.length]  ;
          double[]  dftxn_im=new double[xr.length]  ;
          double[]  dftXk_re=new double[xr.length] ;
          double[]  dftXk_im=new double[xr.length] ;
          double[] idftXK_re=new double[xr.length] ;
          double[] idftXK_im=new double[xr.length] ;
          double[][] idftXN=new double[xr.length][2] ;

     
      if (xr.length==0)
      {
      System.out.println("\n Error");
     
      idftXN[0][0] =0;
      idftXN[0][1]=0;
      return idftXN;
      }
      //xnreal=real(xr);
     
     
      for (int i=0;i<xr.length;i++)
      {
       dftxn_re[i]=xr[i];
       dftxn_im[i]=0;
      }

     //fft
     // dftXk=dft(dftxn,n);

     for (int k=0;k<xr.length;k++)
      {
      dftXk_re[k]=0;
      dftXk_im[k]=0;
      for(int n=0;n<xr.length;n++)
        {
        angle=-2*Math.PI*k*n/xr.length;
        dftXk_re[k]+= dftxn_re[n]*Math.cos(angle)- dftxn_im[n]*Math.sin(angle);
        dftXk_im[k]+= dftxn_re[n]*Math.sin(angle)+ dftxn_im[n]*Math.cos(angle);
        }

     }

     //X={Xk.re,Xk.im};
     
     for (int i=0;i<=xr.length;i++)
      {
      h[i]=0;
      }
    //int no2=(int) Math.floor (xr.length/2);

      if((xr.length>0)&&((2* Math.floor(xr.length/2))==(xr.length))) // n is even
      {

      h[0]=1;
       h[xr.length/2]=1;
       for(int k=1;k<xr.length/2;k++)
       {
       h[k]=2;
       }
      }
      else //n is odd
      {
        if(xr.length>0)
        {
         h[0]=1;
         for(int k=1;k<(xr.length+1)/2;k++)
         {
         h[k]=2;
         }
        }
      }

     for(int i=0;i<xr.length;i++)
      {
      idftXK_re[i]= dftXk_re[i]*h[i];
      idftXK_im[i]= dftXk_im[i]*h[i];
      }


     //idftXN= idft(idftXK, n);
     
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (int n=0;n<xr.length;n++)
      {
      idftXN[n][0] =0;
      idftXN[n][1] =0;
      for(int k=0;k<xr.length;k++)
       {
      
       angle=2*Math.PI*k*n/xr.length;
       idftXN[n][0] +=idftXK_re[k]*Math.cos(angle)-idftXK_im[k]*Math.sin(angle);
       idftXN[n][1] +=idftXK_re[k]*Math.sin(angle)+idftXK_im[k]*Math.cos(angle);
       }
     idftXN[n][0] =idftXN[n][0] /xr.length;
     idftXN[n][1] =idftXN[n][1] /xr.length;

     }


      return idftXN;
     }

       public static  double[][] getanalyticsignal(double [] xn )
         {
        
      int i;
        
            double[][] Hilb =new double[xn.length][2];
            double[][] x_n  =new double[xn.length][2];

     Hilb=hilbert(xn );
     
      for (i=0; i<xn.length; i++)
      {
      x_n[i][0] =xn[i];
       x_n[i][1] =Hilb[i][1] ;
      }
     
     return x_n;
     }

     
      public static  double Arg(double re, double im)
           {
     return Math.atan2(im,re)  ; 
     }

     public static double[][]  get_mod_fi(double[][] x, int n  )
       {
       
       double[][] res=new double[ n][2] ;
       double fi;
       for (int i=0;i<n;i++)
       {
       res[i][0] =Math.sqrt(x[i][0] *x[i][0] +x[i][1] *x[i][1] );
       res[i][1] =Arg(x[i][0] ,x[i][1] );
         { res[i][1] =res[i][1] *180/Math.PI;  }        
       }
       return  res;
     }


     public static double[][]  GetSpectrum(   double [] xn  )
     {

      
      
             double[][] xn_iq=new double[xn.length][2];
      double[][] Xf_IQ =new double[xn.length][2] ;
      double[][] Xf_modfi =new double[xn.length][2]; 
     
      xn_iq =getanalyticsignal( xn );
      Xf_IQ =cdft(xn_iq, xn.length);
      Xf_IQ =ScaleDFT(Xf_IQ,xn.length);
      Xf_modfi =get_mod_fi(Xf_IQ, xn.length );
     
     
            
              
     return   Xf_modfi ;
     }


     public static double[][][]  GetSpectrumarray(   double [] xn , int num,int nseq)
    {
      //int np=(int) num/nseq;
      int npages=(int) Math.floor (num/nseq);

      double [] xnseq =new double [nseq];
      
      double [][]  Xf = new double  [nseq][2];
      double [][][] Xfarray= new double [npages][nseq][2];


     for (int i=0;i<npages; i++) //0,1,32
     {
      for (int j=0;j<nseq;j++)
       {
       xnseq[j]=xn[i+j];     //if  nseq=16, xnseq[j]=xn[k],k= 0...15,16...31
       // if ((i+j)<num) &&( (nseq*npages)>(num/nseq) ) xnseq[j]=0;
       }

      Xf=GetSpectrum(xnseq );
      for (int j=0;j<nseq; j++) //0,1,32
      {
      Xfarray[i][j][0]= Xf[j][0];
      Xfarray[i][j][1]=Xf[j][1];
      }

     }


    return Xfarray;
    }

    public static double get_max(double a,double b )
    {
    if (a>b) {return a;}

    return b;
    }


     public static double[] GetPeakSpectrumComponents(double [] xn , int num, int nseq)
    {
     //int np=(int)num/nseq;
     int npages=(int) Math.floor (num/nseq);
     double [][][] Xfarray =new double [npages][nseq][2];
     double [] Xfpeak=new double[nseq];
     Xfarray=GetSpectrumarray( xn , num, nseq);

     for (int j=0;j<nseq;j++) // for each Xn [frequency](page)
     {
     Xfpeak[j]=0;
     for (int i=0;i<npages;i++) //for pages of spectrum
       {
       Xfpeak[j]=get_max(Xfpeak[j], Xfarray[i][j][0] );
       }

     }

    return Xfpeak;

    }

      public static void main(String[] args) throws IOException
     {

                        Scanner sc = new Scanner(System.in);

       
                     int NMAX=32767;
                     int N=0,n,i, ns=0;
                    double[] xn =new double[NMAX] ;
                    
                      double fs=0,T;
                             
             System.out.print ("\n For input  x[n] from file input '1' ");
      System.out.print ("\n For manual  input  x[n] input   '0' ");
      int key=sc.nextInt();
             
      switch( key)
            {

     

      case 1:
     
               System.out.print ("\n  Open xinput.txt \n");
               String fileNameDefined = "xinputns.txt";
               File file1 = new File(fileNameDefined);
               
            try{
                // -read from filePooped with Scanner class
                Scanner inputStream  = new Scanner(file1);
                fs = (double) inputStream.nextDouble();
                N =inputStream.nextInt();
                ns =inputStream.nextInt(); 
                 System.out.print ("\n ns=");
                 System.out.print (ns);
               for (n=0 ; n<N ;n++)
        {
                 inputStream.nextInt();
              
        xn[n]= inputStream.nextDouble();

        }

               
                
                inputStream.close();


            } catch (FileNotFoundException e) {

                e.printStackTrace();
            }


      break;

      case 0:

       System.out.print ("\n Input fs input,Hz  ");
       fs = (double) sc.nextDouble();
       System.out.print ("\n Input number of samples of xin[n]  ");
          N =sc.nextInt();
                System.out.print ("\n Input number of samples in CDFT sequence  ");
          ns=sc.nextInt();
                 xn =new double[N] ;
       for (n=0 ; n<N ;n++)
       {
       T=n*1/(N *fs );
       System.out.print ("\n n=");
       System.out.print (n);
       System.out.print (" t=");
       System.out.print (T);
       System.out.print (" s; xn[n]=  " ); 
       xn[n]= sc.nextDouble();
       }
        break;
      default: 
          break;
     }
      System.out.print ("\n  fs =");
       System.out.print (fs);
      System.out.print (" Hz");
       System.out.print ("\n Nsamples1=");
       System.out.print (N);

     for (n=0 ; n<N ;n++)
     {
      T=n*1/(N *fs );
       System.out.print ("\n n=");
       System.out.print (n);
       System.out.print (" x[n]=");
       System.out.print (xn[n]);
       System.out.print ("  t= ");
       System.out.print (T);
       System.out.print (" s");
     }
     
       System.out.println ("\n Input 0 for continue ");
      sc.nextLine();
       System.out.print ("\n OK");

     

                      System.out.println("Array of Xn is: ");
        for (  i=0;i<N;i++)
                      {
                       System.out.print("\n");
                       System.out.print("N= ");
                       System.out.print(i);                
                       System.out.print (" ; xn[n]= ");
                       System.out.print(xn[i]);
                      }

                     double Xfmax[]  =new  double[ns]  ;
                    System.out.println("\n Parsing spectrum ");
                    Xfmax=GetPeakSpectrumComponents(xn, N, ns );
                     //Xf=GetSpectrum(xn, N);


                  
                System.out.println("\n Array of Xfmax is: ");
        for (i=0;i<ns ;i++)
                      {
                       System.out.print("\n");
                       System.out.print("ns= ");
                       System.out.print(i);                

                        System.out.print (" ; max |Xf[k]|= ");
                        System.out.print(Xfmax[i]   );
                        
                        System.out.print (" ; f=");
                        System.out.print ( fs*i/ns );
                        System.out.print(" Hz");

                      }

        File outFile = new File ("maxspectr.txt");
        FileWriter fWriter = new FileWriter (outFile);
        PrintWriter pWriter = new PrintWriter (fWriter);


        pWriter.println("\n");
        pWriter.print(String.format("%e",fs ));
        pWriter.print("   ");
        pWriter.print(String.format("%d",N ) );
                       
        pWriter.println("\n");


       for (i=0; i<ns ;i++)
       { 

     
             pWriter.println("\n ");
      pWriter.print(String.format("%d",i) );
             pWriter.print("   "); 
      pWriter.print(String.format("%e",Xfmax[i] ));
      pWriter.print("  ");
             pWriter.print (String.format("%e", fs*i/ns) );
             pWriter.print("  ");    
       }
     
     
        pWriter.close();

                    System.out.println("\n End of sequence ");

                       sc.close();
     }

    }

     

    Friday, January 26, 2018 3:01 PM
  • You may convert Java programs into VC++ console applications (if the problem of work with big 2-dimension arrays is solved ).
    Friday, January 26, 2018 3:04 PM
  • Program for obtaining envelope of sinal sequence using Hilbert transform .

    // getenv.cpp : Defines the entry point for the console application.
    //

    #include "stdafx.h"

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


    #define NMAX 1000
    #define M_PI 3.1415926535897932384626433832795

    typedef struct {
    double re[NMAX];
    double im[NMAX];
    } CData;

     typedef struct  {
         double data[NMAX];
     } Array;


    CData  cdft(CData xn  ,int N )
     {
       
      
     CData Xk ;
     int k,n  ;
     double angle;

     for (k=0;k<N;k++)
      {
      Xk.re[k]=0;
      Xk.im[k]=0;
      for(n=0;n<N;n++)
        {
        angle=-2*M_PI*k*n/N;  
        Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
        Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
        }

     }

     return Xk;
     }

     CData  cidft(CData idft_Xk, int N  )
     {

     int k,n ;
     double angle;
     CData  idft_xn;
     
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (n=0;n<N;n++)
      {
      idft_xn.re[n]=0;
      idft_xn.im[n]=0;
      for(k=0;k<N;k++)
       {
     
       angle=2*M_PI*k*n/N;
       idft_xn.re[n]=idft_xn.re[n]+idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
       idft_xn.im[n]=idft_xn.im[n]+idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
       }

     idft_xn.im[n]=idft_xn.im[n]/N;
     idft_xn.re[n]=idft_xn.re[n]/N;

     }

    return idft_xn ;
     }

    CData ScaleDFT(CData xk,int  num)
     {

     int i;
     for (i=0; i<num ; i++)
       {
        xk.re[i]=xk.re[i]/num;
        xk.im[i]=xk.im[i]/num;
       }
     return xk;
     }

    CData ScaleIQidft(CData Xk, int num)
     {
     int i;
       for (i=0;  i<num ; i++)
      {
       Xk.re[i]=Xk.re[i]*num;
       Xk.im[i]=Xk.im[i]*num;
      }

    return Xk;

    }

     CData hilbert(double *xr , int N)
     {
      int  k,n,i, no2 ;
      double angle;
      double h[NMAX];
      double dft__xn_re[NMAX];
      double dft__xn_im[NMAX];
       double dft__Xk_re[NMAX];
       double dft__Xk_im[NMAX];
         double idft__XK_re[NMAX];
         double idft__XK_im[NMAX];
       CData idft__XN;

     
      if (N==0)
      {
      printf("\n Error");
      idft__XN.re[0]=0;
      idft__XN.im[0]=0;
      return idft__XN;
      }
      //xnreal=real(xr);
      no2=(int)floor ((double)N/2);
     
      for (i=0;i<N;i++)
      {
       dft__xn_re[i]=xr[i];
       dft__xn_im[i]=0;
      }

     //fft
     // dft__Xk=dft(dft__xn,n);

    for (k=0;k<N;k++)
      {
      dft__Xk_re[k]=0;
      dft__Xk_im[k]=0;
      for(n=0;n<N;n++)
        {
        angle=-2*M_PI*k*n/N;
        dft__Xk_re[k]+= dft__xn_re[n]*cos(angle)- dft__xn_im[n]*sin(angle);
        dft__Xk_im[k]+= dft__xn_re[n]*sin(angle)+ dft__xn_im[n]*cos(angle);
        }

     }

     //X={Xk.re,Xk.im};

     for (i=0;i<=N;i++)
      {
      h[i]=0;
      }


      if((N>0)&&((2*floor((double)N/2))==(N))) // n is even
      {

      h[0]=1;
       h[N/2]=1;
       for(k=1;k<N/2;k++)
       {
       h[k]=2;
       }
      }
      else //n is odd
      {
        if(N>0)
        {
         h[0]=1;
         for(k=1;k<(N+1)/2;k++)
         {
         h[k]=2;
         }
        }
      }

     for(i=0;i<N;i++)
      {
      idft__XK_re[i]= dft__Xk_re[i]*h[i];
      idft__XK_im[i]= dft__Xk_im[i]*h[i];
      }


     //idft__XN= idft(idft__XK, n);
     
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (n=0;n<N;n++)
      {
      idft__XN.re[n]=0;
      idft__XN.im[n]=0;
      for(k=0;k<N;k++)
       {
       //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
       angle=2*M_PI*k*n/N;
       idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
       idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
       }
     idft__XN.im[n]=idft__XN.im[n]/N;
     idft__XN.re[n]=idft__XN.re[n]/N;

     }


    return idft__XN;
     }

    CData getanalyticsignal(double *xn, int num)
     {
     CData x_n,Hilb;
     int i;
       //x=x+jy,y=0
       //hilb(x)=a+jb
       //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
       //I=x(t), Q=j*Hilb(xre(t))

    Hilb=hilbert(xn,num);
     
     for (i=0; i<num; i++)
     {
      x_n.re[i]=xn[i];
      x_n.im[i]=Hilb.im[i];
     }
     
     return x_n;
     }


     
    double Arg(double re, double im)
    {
     
    return atan2 (im,re)  ; 
    }

    CData  get_mod_fi(CData x, int n, int mode   )
     {
      int i;
      CData res;
      double fi;
      for (i=0;i<n;i++)
      {
      res.re[i]=pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
      //double treshh=1e9;
      // x.im[i]=x.im[i]*treshh;
      // x.im[i]=floor(x.im[i])/treshh;
      // x.re[i]=x.re[i]*treshh;
      // x.re[i]=floor(x.re[i])/treshh;
      res.im[i]=Arg(x.re[i],x.im[i]);
      if (mode==1)  { res.im[i]=res.im[i]*180/M_PI;  }
              
      }
      return  res;
     }

    CData GetSpectrum(double *xn, int num)
    {
     
    CData xn_iq, Xf_IQ, Xf_modfi;
    xn_iq=getanalyticsignal( xn,  num);
     
    Xf_IQ=cdft(xn_iq  ,num );
    Xf_IQ=ScaleDFT(Xf_IQ,num);
    Xf_modfi=get_mod_fi(Xf_IQ, num,1   );
     
    return Xf_modfi;
    }

    Array getenvelope(double *x, int num)
     {
     int i;
     CData  y;
     Array env;
     y=hilbert(x,num);
     for (i=0;i<num;i++)
      {
      env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
      }
      return env;
     }
     
     
     
    double getmax( Array   k , int N)
    {
    double kmax=k.data[0];
    double k1=0;
    int i;
    for (i=1;i<N;i++)
    {
     k1=k.data[i];
      if(fabs( k1 )>fabs(kmax)) { kmax= k.data[i]; } 
      
    }
     
     return kmax;
    }
     
    double getmin( Array  k, int N)
    {
    double kmin=k.data[0];
    int i;
    for (i=0;i<N;i++)
    {
      if(fabs(k.data[i])<fabs(kmin)) { kmin= k.data[i]; } 
    }
     return kmin;
    }

    double getmid( Array  k,int N)
    {
    double kmid=0;
    int i;
    for (i=0;i<N;i++)
    {
     kmid+=   k.data[i];
    }
    kmid=kmid/N;
     return kmid;
    }


    FILE   *fp1,*fp2 ;

    int _tmain(int argc, _TCHAR* argv[])
    //int main ()
    {
    double delta =1e-15;
    double xin[NMAX];
    double yout[NMAX];
    Array  Xenv1  ;
    double Kmax=0, Kmin=0, Kmid=0;
    int key;
    int k,n,i,N1,N2;
    double fs1,fs2,T,x,f;


    /***********input x[n] *************/
    label_00:
    printf ("\n For input  x[n] from file input '1' ");
    printf ("\n For manual  input  x[n] input   '0' ");
    std::cin>>key;
    if (key== 1  ) { goto label_01; }
    if (key== 0 ) { goto label_02; }

    label_01:

     printf("\n Open xinput.txt for reading x[n]   ");
    if ((fp1=fopen("xinput.txt","r"))==NULL)  printf ("\n File opening error ");
    fscanf(fp1,"%le",&fs1);
    fscanf(fp1,"%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     fscanf (fp1,"%d",&k);
     fscanf (fp1,"%le",&xin[n]);
    }


    fclose(fp1);

    goto label_03;

    label_02:

    printf("\n Input fs input,Hz  " ) ;
    scanf( "%le",&fs1);
    printf("\n Input number of samples of xin[n]  ") ;
    scanf( "%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     T=n*1/(N1*fs1); 
     printf ("\n n=%d ; t=%le s; xin[n]= ",n,T);
     scanf ( "%le",&xin[n]);
    }

     label_03:
    printf( "\n fs1=%le Hz ",fs1);
    printf( "\n Nsamples1= %d ",N1);

    for (n=0 ; n<N1;n++)
    {
     T=n*1/(N1*fs1); 
     printf("\n n=%d x[n]=%le,  t=%le s",n,xin[n],T);
    }

    printf("\n Input 0 for continue ");
    std::cin>>key;
    printf("\n OK.   ");

     


    /*****************************/

    label1:
     
    printf("\n Obtaining envelope of x[n]");
    Xenv1=getenvelope(xin, N1);

    printf("\n |Hilbert(x[n])| : ");
     
    if ((fp2=fopen("envelope.txt","a"))==NULL)  printf ("\n File opening error");
    fprintf(fp2,"\n\n   %le   %d \n",fs1, N1);


    for (i=0; i<N1;i++)
    {
     
     
    printf ("\n %d   |Hilbert(x[n])|= %le  ",  i, Xenv1.data[i] );
    fprintf (fp2,"\n %d   %le  ",  i, Xenv1.data[i]);

    }
    fclose(fp2);

    printf ("\n  input   '0' for continue ");
    std::cin>>key;

    Kmax=getmax(Xenv1,N1);
    Kmin=getmin(Xenv1,N1);
    Kmid=getmid(Xenv1,N1);
    printf ("\n  Kmax= %le ",   Kmax);
    printf ("\n  Kmin= %le ",   Kmin);
    printf ("\n  Kmid= %le ",   Kmid);

    if ((fp2=fopen("coeffs.txt","a"))==NULL)  printf ("\n File opening error");
    fprintf(fp2,"\n\n   %le   %d \n",fs1,N1);
    fprintf (fp2,"\n  Kmax= %le ",   Kmax);
    fprintf (fp2,"\n  Kmin= %le ",   Kmin);
    fprintf (fp2,"\n  Kmid= %le ",   Kmid);
    fclose(fp2);

    printf ("\n  input   '0' for exit ");
    std::cin>>key;

    return 0;

    }

    // stdafx.cpp : source file that includes just the standard includes
    // getenv.pch will be the pre-compiled header
    // stdafx.obj will contain the pre-compiled type information

    #include "stdafx.h"

    // TODO: reference any additional headers you need in STDAFX.H
    // and not in this file

    // stdafx.h : include file for standard system include files,
    // or project specific include files that are used frequently, but
    // are changed infrequently
    //

    #pragma once


    #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
    #include <stdio.h>
    #include <tchar.h>

    // TODO: reference additional headers your program requires here

    Friday, January 26, 2018 4:19 PM
  • Program for obtaining of the Bessel's  functions for obtainer of the PWM and FM spectrums of  signals.

    // bessel.cpp : Defines the entry point for the console application.
    //

    #include "stdafx.h"


    #include <stdio.h>
    #include <stdlib.h>
    #include <iostream>
    #include <math.h>
    #include "bessel.h"


    int main()

    {
    int n;
    double x;
    double result1, result2,result3,result4;

    std::cout<<"Input order n : "<<"\n";
    std::cin>>n;
    std::cout<<"Input x : "<<"\n";
    std::cin>>x;
    result1=jn(n,x);
    std::cout<<"jn(n,x)= "<<result1<<"\n";
    result2=bessj( n, x );
    std::cout<<"jn(n,x)= "<<result2<<"\n";


    result3=yn(n,x);
    std::cout<<"yn(n,x)= "<<result3<<"\n";
    result4=bessy(  n,  x );
    std::cout<<"yn(n,x)= "<<result4<<"\n";

    std::cout<<"Input '0' for exit";
    char key;
    std::cin>>key;
    return 0;
    }

    //bessel.h (add into project using Add Item-> .h file,then  put the code into file)

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


    #define ACC 40.0
    #define BIGNO 1.0e10
    #define BIGNI 1.0e-10

    /*
    http://www.astro.rug.nl/~gipsy/sub/bessel.c
    */

    static double bessj0( double x )
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate Bessel function of first kind and order  */
    /*          0 at input x                                      */
    /*------------------------------------------------------------*/
    {
       double ax,z;
       double xx,y,ans,ans1,ans2;

       if ((ax=fabs(x)) < 8.0) {
          y=x*x;
          ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
             +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
          ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
             +y*(59272.64853+y*(267.8532712+y*1.0))));
          ans=ans1/ans2;
       } else {
          z=8.0/ax;
          y=z*z;
          xx=ax-0.785398164;
          ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
             +y*(-0.2073370639e-5+y*0.2093887211e-6)));
          ans2 = -0.1562499995e-1+y*(0.1430488765e-3
             +y*(-0.6911147651e-5+y*(0.7621095161e-6
             -y*0.934935152e-7)));
          ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
       }
       return ans;
    }

    static double bessj1( double x )
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate Bessel function of first kind and order  */
    /*          1 at input x                                      */
    /*------------------------------------------------------------*/
    {
       double ax,z;
       double xx,y,ans,ans1,ans2;

       if ((ax=fabs(x)) < 8.0) {
          y=x*x;
          ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
             +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
          ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
             +y*(99447.43394+y*(376.9991397+y*1.0))));
          ans=ans1/ans2;
       } else {
          z=8.0/ax;
          y=z*z;
          xx=ax-2.356194491;
          ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
             +y*(0.2457520174e-5+y*(-0.240337019e-6))));
          ans2=0.04687499995+y*(-0.2002690873e-3
             +y*(0.8449199096e-5+y*(-0.88228987e-6
             +y*0.105787412e-6)));
          ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
          if (x < 0.0) ans = -ans;
       }
       return ans;
    }

    double bessj( int n, double x )
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate Bessel function of first kind and order  */
    /*          n at input x                                      */
    /* The function can also be called for n = 0 and n = 1.       */
    /*------------------------------------------------------------*/
    {
       int    j, jsum, m;
       double ax, bj, bjm, bjp, sum, tox, ans;


       if (n < 0)
       {
          double   dblank;
          printf("Error");
        //  setdblank_c( &dblank );
          return( 0 );
       }
       ax=fabs(x);
       if (n == 0)
          return( bessj0(ax) );
       if (n == 1)
          return( bessj1(ax) );

       if (ax == 0.0)
          return 0.0;
       else if (ax > (double) n) {
          tox=2.0/ax;
          bjm=bessj0(ax);
          bj=bessj1(ax);
          for (j=1;j<n;j++) {
             bjp=j*tox*bj-bjm;
             bjm=bj;
             bj=bjp;
          }
          ans=bj;
       } else {
          tox=2.0/ax;
          m=2*((n+(int) sqrt(ACC*n))/2);
          jsum=0;
          bjp=ans=sum=0.0;
          bj=1.0;
          for (j=m;j>0;j--) {
             bjm=j*tox*bj-bjp;
             bjp=bj;
             bj=bjm;
             if (fabs(bj) > BIGNO) {
                bj *= BIGNI;
                bjp *= BIGNI;
                ans *= BIGNI;
                sum *= BIGNI;
             }
             if (jsum) sum += bj;
             jsum=!jsum;
             if (j == n) ans=bjp;
          }
          sum=2.0*sum-bj;
          ans /= sum;
       }
       return  x < 0.0 && n%2 == 1 ? -ans : ans;
    }


    static double bessy0( double x )
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate Bessel function of second kind and order */
    /*          0 at input x.                                     */
    /*------------------------------------------------------------*/
    {
       double z;
       double xx,y,ans,ans1,ans2;

       if (x < 8.0) {
          y=x*x;
          ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
             +y*(10879881.29+y*(-86327.92757+y*228.4622733))));
          ans2=40076544269.0+y*(745249964.8+y*(7189466.438
             +y*(47447.26470+y*(226.1030244+y*1.0))));
          ans=(ans1/ans2)+0.636619772*bessj0(x)*log(x);
       } else {
          z=8.0/x;
          y=z*z;
          xx=x-0.785398164;
          ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
             +y*(-0.2073370639e-5+y*0.2093887211e-6)));
          ans2 = -0.1562499995e-1+y*(0.1430488765e-3
             +y*(-0.6911147651e-5+y*(0.7621095161e-6
             +y*(-0.934945152e-7))));
          ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
       }
       return ans;
    }

    static double bessy1( double x )
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate Bessel function of second kind and order */
    /*          1 at input x.                                     */
    /*------------------------------------------------------------*/
    {
       double z;
       double xx,y,ans,ans1,ans2;

       if (x < 8.0) {
          y=x*x;
          ans1=x*(-0.4900604943e13+y*(0.1275274390e13
             +y*(-0.5153438139e11+y*(0.7349264551e9
             +y*(-0.4237922726e7+y*0.8511937935e4)))));
          ans2=0.2499580570e14+y*(0.4244419664e12
             +y*(0.3733650367e10+y*(0.2245904002e8
             +y*(0.1020426050e6+y*(0.3549632885e3+y)))));
          ans=(ans1/ans2)+0.636619772*(bessj1(x)*log(x)-1.0/x);
       } else {
          z=8.0/x;
          y=z*z;
          xx=x-2.356194491;
          ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
             +y*(0.2457520174e-5+y*(-0.240337019e-6))));
          ans2=0.04687499995+y*(-0.2002690873e-3
             +y*(0.8449199096e-5+y*(-0.88228987e-6
             +y*0.105787412e-6)));
          ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
       }
       return ans;
    }


    double bessy( int n, double x )
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate Bessel function of second kind and order */
    /*          n for input x. (n >= 0)                           */
    /* Note that for x == 0 the functions bessy and bessk are not */
    /* defined and a blank is returned.                           */
    /*------------------------------------------------------------*/
    {
       int j;
       double by,bym,byp,tox;


       if (n < 0 || x == 0.0)
       {
          double   dblank;
        printf("Error"); // setdblank_c( &dblank );
          return(0 );
       }
       if (n == 0)
          return( bessy0(x) );
       if (n == 1)
          return( bessy1(x) );

       tox=2.0/x;
       by=bessy1(x);
       bym=bessy0(x);
       for (j=1;j<n;j++) {
          byp=j*tox*by-bym;
          bym=by;
          by=byp;
       }
       return by;
    }

    static double bessi0( double x )
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate modified Bessel function In(x) and n=0.  */
    /*------------------------------------------------------------*/
    {
       double ax,ans;
       double y;


       if ((ax=fabs(x)) < 3.75) {
          y=x/3.75,y=y*y;
          ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
             +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
       } else {
          y=3.75/ax;
          ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
             +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
             +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
             +y*0.392377e-2))))))));
       }
       return ans;
    }


    static double bessi1( double x)
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate modified Bessel function In(x) and n=1.  */
    /*------------------------------------------------------------*/
    {
       double ax,ans;
       double y;


       if ((ax=fabs(x)) < 3.75) {
          y=x/3.75,y=y*y;
          ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
             +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
       } else {
          y=3.75/ax;
          ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
             -y*0.420059e-2));
          ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
             +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
          ans *= (exp(ax)/sqrt(ax));
       }
       return x < 0.0 ? -ans : ans;
    }


    double bessi( int n, double x)
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate modified Bessel function In(x) for n >= 0*/
    /*------------------------------------------------------------*/
    {
       int j;
       double bi,bim,bip,tox,ans;


       if (n < 0)
       {
          double   dblank;
         printf("Error");// setdblank_c( &dblank );
          return( 0);
       }
       if (n == 0)
          return( bessi0(x) );
       if (n == 1)
          return( bessi1(x) );


       if (x == 0.0)
          return 0.0;
       else {
          tox=2.0/fabs(x);
          bip=ans=0.0;
          bi=1.0;
          for (j=2*(n+(int) sqrt(ACC*n));j>0;j--) {
             bim=bip+j*tox*bi;
             bip=bi;
             bi=bim;
             if (fabs(bi) > BIGNO) {
                ans *= BIGNI;
                bi *= BIGNI;
                bip *= BIGNI;
             }
             if (j == n) ans=bip;
          }
          ans *= bessi0(x)/bi;
          return  x < 0.0 && n%2 == 1 ? -ans : ans;
       }
    }


    static double bessk0( double x )
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate modified Bessel function Kn(x) and n=0.  */
    /*------------------------------------------------------------*/
    {
       double y,ans;

       if (x <= 2.0) {
          y=x*x/4.0;
          ans=(-log(x/2.0)*bessi0(x))+(-0.57721566+y*(0.42278420
             +y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2
             +y*(0.10750e-3+y*0.74e-5))))));
       } else {
          y=2.0/x;
          ans=(exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1
             +y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
             +y*(-0.251540e-2+y*0.53208e-3))))));
       }
       return ans;
    }


    static double bessk1( double x )
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate modified Bessel function Kn(x) and n=1.  */
    /*------------------------------------------------------------*/
    {
       double y,ans;

       if (x <= 2.0) {
          y=x*x/4.0;
          ans=(log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144
             +y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1
             +y*(-0.110404e-2+y*(-0.4686e-4)))))));
       } else {
          y=2.0/x;
          ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619
             +y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2
             +y*(0.325614e-2+y*(-0.68245e-3)))))));
       }
       return ans;
    }

    double bessk( int n, double x )
    /*------------------------------------------------------------*/
    /* PURPOSE: Evaluate modified Bessel function Kn(x) and n >= 0*/
    /* Note that for x == 0 the functions bessy and bessk are not */
    /* defined and a blank is returned.                           */
    /*------------------------------------------------------------*/
    {
       int j;
       double bk,bkm,bkp,tox;


       if (n < 0 || x == 0.0)
       {
          double   dblank;
         printf("Error");// setdblank_c( &dblank );
          return( 0 );
       }
       if (n == 0)
          return( bessk0(x) );
       if (n == 1)
          return( bessk1(x) );

       tox=2.0/x;
       bkm=bessk0(x);
       bk=bessk1(x);
       for (j=1;j<n;j++) {
          bkp=bkm+j*tox*bk;
          bkm=bk;
          bk=bkp;
       }
       return bk;
    }


    #undef ACC
    #undef BIGNO
    #undef BIGNI

    // stdafx.cpp : source file that includes just the standard includes
    // bessel.pch will be the pre-compiled header
    // stdafx.obj will contain the pre-compiled type information

    #include "stdafx.h"

    // TODO: reference any additional headers you need in STDAFX.H
    // and not in this file

    // stdafx.h : include file for standard system include files,
    // or project specific include files that are used frequently, but
    // are changed infrequently
    //

    #pragma once


    #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
    #include <stdio.h>
    #include <tchar.h>
    //#include "bessel.h"


    // TODO: reference additional headers your program requires here


    • Edited by USERPC01 Friday, January 26, 2018 6:57 PM bug fix
    Friday, January 26, 2018 6:54 PM
  • project for obtaining IR coefficients of FIR from AFR data   with #include <vector>,  #include <complex>, #include <cmath>

    //afrtoht.cpp

    // afrtoht.cpp : main project file.

    #include "stdafx.h"
    #include <stdio.h>
    /*
    #include <tchar.h>
     #include <iostream>
    #include <stdlib.h>
    #include <cmath>
    */
    #include <iostream>
    #include <fstream>
    #include <complex>
    #include <vector>
    #include <set>
    #include <cmath>
    #include <limits>

    using namespace std;

    #define M_PI 3.1415926535897932384626433832795
    //int _tmain(int argc, _TCHAR* argv[])
    //{
    // return 0;
    //}

     //afrtohtcoefs.cpp
    /* Define complex double data type */
     typedef complex<double> dcomp;

     //std::vector<std::vector<dcomp> >
     typedef std::vector<dcomp> cplxvector ;

    #define  NMAX 512
    //FILE *fp, *fp1;
     
    int k,n;
    int w,t,N;
    double h[NMAX];
    double h_im[NMAX];
    double Re_H[NMAX];
    double Im_H[NMAX];
    double Re_x[NMAX];
    double Im_x[NMAX];
    double abs_H[NMAX];
    double arg_H[NMAX];
    double f[NMAX];
    float angle;
    double   K, fs,T,finp,Tmax;
    char str;
    char key;

    /*
    float heavi( double t)
    {
    if (t<0 ) return 0;
    if (t>=0) return 1;
    }

    float heavi1( int n)
    {
    if (n<0 ) return 0;
    if (n>=0) return 1;
    }*/

    double getarg(float Re_H,float Im_H)
    {

    return atan2(Im_H,Re_H) * 180 / M_PI;
    }

    cplxvector CIDFT(cplxvector idftXk, int N)
    {

     
     double angle;
     cplxvector  idftxn;
     dcomp data;

    //N=sizeof(xk)/sizof(double);
     //ifft
      for (int n=0;n<N;n++)
      {
      data=complex<double>(0.0,0.0);
      
     
      for(int k=0;k<N;k++)
       {
        angle=2*M_PI*k*n/N;
        data+=idftXk.at(k) *  exp(complex<double> (0.0,angle));
       // data+=idft_Xk.at(k)* cmath.exp( (0.0,angle));
       //data+=idft_Xk.at(k)*  exp( 1i*angle );
       //idft_xn.re[n]+= idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
       // idft_xn.im[n]+= idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
        }

    data/=N;
    idftxn.push_back( data);

     }
    return idftxn;
    }


     int   main( )
    {

     fstream finp;
     std::cout <<  "\n For input  |H(jw)|, fi(w) from file input '1' " ;
     std::cout <<  "\n For manual  input |H(jw)|, fi(w) input   '0' " ;

        std::cin>>key;
      ofstream fout ("hcoefs1.txt" , ios::ate | ios::app);
         fout.close();
       if (key=='0' ) { goto label_02;}
       if (key=='1' ) { goto label_01;}

    label_01:;
     std::cout << "\n Open afrcoefs.txt for reading AFR coefficients ";
     
    // ifstream finp("afrcoefs.txt", ios::in );
     finp.open("afrcoefs.txt", ios::in );
     
      finp>>fs ;
      finp>>N ;
    Tmax=N*1/fs;
    for (k=0 ; k<N;k++)
    {
    // f[k]=k*fs/N;
     
      //fout<<"\n n="<<k<<" |H(jw)|="<<abs_H[k]<< "; arg(H(jw))= "<< arg_h[k]<<"; f="<<f[k]<<" Hz"  ;
      finp>> n ;
      finp>>abs_H[k] ;
      finp>>arg_H[k] ;
     
    //  finp>>f[k] ;
     
     arg_H[k]*= M_PI/180;
    }
    finp.close();
    fout.close();
     
    goto label_03;


    label_02:;

    std::cout<<"\n Input fs ,Hz " ;
    std::cin>>fs;
    std::cout<<" Input  number of points " ;
    std::cin>>N ;
    Tmax=N*1/fs;
    std::cout<<"\n Tmax ="<<Tmax<<" s" ;
    for (k=0 ; k<N;k++)
    {
     f[k]=k*fs/N;

    //fout<<"\n n="<<k<<"  |H(jw)|=" <<  abs_H[k]<<" ; arg(H(jw))="<< arg_h[k]<<"; f="<<f[k]<<" Hz" ;
    //fout<<"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
       std::cout<<"\n F="<<f[k]<<" Hz" ;
       std::cout<<"\n k="<<k<<";    |H(jw)|[k]=" ;
       std::cin>> abs_H[k] ;
       std::cout<<" k="<<k<<" ;   arg(H(jw))[k]=" ;
        std::cin>>arg_H[k] ;
      //abs_H[k]=1;
      //arg_H[k]=0;
    arg_H[k]*= M_PI/180;

    }


    label_03:
     
    for (k=0 ; k<N;k++)
    {
    f[k]=k*fs/N;
     std::cout << "\n n="<<k<<";|H(jw)|= "<<abs_H[k]<< " ;arg(H(jw))="<<arg_H[k]*180/M_PI<<"  deg;f="<< f[k]<<" Hz";
     
    }


    std::cout<< "\n Input 0 for continue " ;
    std::cin.get();
    std::cout<<"\n OK. Parsing ..." ;
     
    // H(jw)=integr(0;inf;h(tau)*exp(-j*w*tau) ; d tau)
    //K=2*3.1415926/fs;
     /*
    DFT
    X[k]=sum(n=0;N-1; x[n]*exp(-2*pi*j*k*n/N) ), k=0,...,N-1
     IDFT
    x[n]=(1/N)*sum(k=0;N-1;X[k]*exp(2*pi*k*n/N) ), n=0,...,N-1
     100*10^-6*500
    */

    //label1:;


    //cidft
    //printf("\n Open hcoefs1.txt for writing coeffs h[n] ");


     
     
     cplxvector  idft_xn, idft_Xk;
     dcomp data;

     for(int k=0;k<N;k++)
      {
       Re_H[k]=abs_H[k]*cos(arg_H[k]);
       Im_H[k]=abs_H[k]*sin(arg_H[k]);
       data =complex<double> (Re_H[k] ,Im_H[k]) ; 
      // std::cout<<"\n k="<<k<< "  Xk[k]=  " <<data  ; /*real(idft_xn.at(n) ) */   ;
       idft_Xk.push_back(data); 
     }


     std::cout <<"\n IDFT\n";
     idft_xn= CIDFT(idft_Xk,N);

    fstream fout1 ("hcoefs1.txt" , ios::ate | ios::app);
      for (int n=0;n<N;n++)
      {
      std::cout<<"\n n="<<n<< "  h[n]=  " <<real(idft_xn.at(n) )  ;
      fout1 << "\n"<<n <<"  "<< real(idft_xn.at(n) )   ; 
      }
       fout1.close();

     //fclose (fp1);
     std::cout<< "\n  Re(h(n)=h[n]), use it to load into FIR  "   ;
    //label3:;
    std::cout<<"\n Input 0 for exit " ;
    std::cin.get();
    return 0;
    }

    //stdafx.cpp

    // stdafx.cpp : source file that includes just the standard includes
    // afrtoht.pch will be the pre-compiled header
    // stdafx.obj will contain the pre-compiled type information

    #include "stdafx.h"

    //stdafx.h

    // stdafx.h : include file for standard system include files,
    // or project specific include files that are used frequently, but
    // are changed infrequently
    //

    #pragma once


     #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers


    // TODO: reference additional headers your program requires here

    //afrcoefs.txt

    44100
    4
    0  1  0
    1  1  0
    2  0  0
    3  0  0

    //hcoefs1.txt


    0  0.5
    1  0.25
    2  0
    3  0.25


    • Edited by USERPC01 Tuesday, February 6, 2018 6:22 PM fix bugs
    Tuesday, February 6, 2018 6:20 PM
  • Example of using <vector>, <complex> for  Hilbert transform

    (fix bugs)

    // cplxarray.cpp : main project file.

    #include "stdafx.h"
     
    #include <iostream>
    //#include <tgmath>//for the type generate marcros.

    #include <complex>
    #include <vector>
    #include <set>
    #include <cmath>
    #include <limits>
    /*
    #include <cstdio>
    #include <cstdlib>
    #include <cmath>
    #include <limits>
    #include <iostream>
     */
    using namespace std;

    /* Define complex double data type */
     typedef complex<double> dcomp;

     //std::vector<std::vector<dcomp> >
     typedef std::vector<dcomp> cplxvector ;

    #define NMAX 2048

    #define M_PI 3.1415926535897932384626433832795
    /*
    typedef struct {

    double re[NMAX];
    double im[NMAX];
    } CData;

     typedef struct  {
         double data[NMAX];
     } Array;

    */

    cplxvector cdft( cplxvector xn  ,int N )
     {
        
     
     cplxvector Xk ;
     dcomp data;
     int k,n  ;
     double angle;
     
     for (k=0;k<N;k++)
      {
      //Xk.re[k]=0;
      //Xk.im[k]=0;
     data=(0.0,0.0);
      for(n=0;n<N;n++)
        {
         angle=-2*M_PI*k*n/N;  
         //Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
         //Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
      //data+=xn.at(n)*cmath.exp( (0.0,angle));

      data+=xn.at(n)* exp( (0.0,angle));
        }

     }

     return Xk;
     }

     cplxvector  cidft(cplxvector idft_Xk, int N  )
     {

     int k,n ;
     double angle;
     cplxvector  idft_xn;
     dcomp data; 
     //N=sizeof(xk)/sizof(double);
     //ifft
      for (n=0;n<N;n++)
      {
      data=(0.0,0.0);
      
     
      for(k=0;k<N;k++)
       {
     
        angle=2*M_PI*k*n/N;
     //creal(z1),
     //cimag(z1),
    // data+=idft_Xk.at(k)* cmath.exp( (0.0,angle));
     data+=idft_Xk.at(k)*  exp( (0.0,angle));
       //data+=idft_Xk.at(k)*  exp( 1i*angle );
       //idft_xn.re[n]+= idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
       // idft_xn.im[n]+= idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
        }

    data/=N;
    idft_xn.push_back( data);
    // idft_xn.im[n]=idft_xn.im[n]/N;
    //  idft_xn.re[n]=idft_xn.re[n]/N;

     }

    return idft_xn ;
     }


     cplxvector  ScaleDFT( cplxvector  xk, int  num)
     {

     int i;
     for (i=0; i<num ; i++)
       {
         xk.at(i)/= num;
      //  xk.im[i]=xk.im[i]/num;
       }
     return xk;
     }

     cplxvector  ScaleIQidft( cplxvector  Xk, int num)
     {
     int i;
        for (i=0;  i<num ; i++)
      {
      Xk.at(i)*=num;
       //Xk.re[i]=Xk.re[i]*num;
       //Xk.im[i]=Xk.im[i]*num;
      }

    return Xk;

    }
     
      cplxvector  hilbert(double *xr , int N)
     {
      int  k,  i, no2 ;
      //double angle;
      double h[NMAX];
      cplxvector dft__xn;
      //double dft__xn_re[NMAX];
      //double dft__xn_im[NMAX];
      cplxvector dft__Xk;
       //double dft__Xk_re[NMAX];
       //double dft__Xk_im[NMAX];
      cplxvector idft__XK;
        // double idft__XK_re[NMAX];
        // double idft__XK_im[NMAX];
        cplxvector  idft__XN;

     
      if (N==0)
      {
      printf("\n Error");
       idft__XN.push_back((0,0)) ;
     
      return idft__XN;
      }
      //xnreal=real(xr);
      no2=(int) floor ((double)N/2);
     // dcomp  tmpdata;

      for (i=0;i<N;i++)
      {
       // tmpdata=(xr[i],0) ;
         dft__xn.push_back((xr[i],0))  ;
       //dft__xn_re[i]=xr[i];
          //dft__xn_im[i]=0;
      }

     //fft
     // dft__Xk=dft(dft__xn,n);
    dft__Xk=cdft( dft__xn  , N );

     //X={Xk.re,Xk.im};

     for (i=0;i<=N;i++)
      {
      h[i]=0;
      }


      if((N>0)&&(( 2*floor((double)N/2) )==(N))) // n is even
      {

      h[0]=1;
       h[N/2]=1;
       for(k=1;k<N/2;k++)
       {
       h[k]=2;
        }
      }
      else //n is odd
      {
         if(N>0)
        {
          h[0]=1;
         for(k=1;k<(N+1)/2;k++)
         {
          h[k]=2;
         }
         }
      }

     for(i=0;i<N;i++)
      {
      idft__XK.at(i)*=h[i];
      //idft__XK_re[i]= dft__Xk_re[i]*h[i];
      //idft__XK_im[i]= dft__Xk_im[i]*h[i];
      }


     //idft__XN= idft(idft__XK, n);
     
      //N=sizeof(xk)/sizof(double);
     //ifft
    idft__XN=cidft(idft__XK,  N  );


     /*
      for (n=0;n<N;n++)
      {
      idft__XN.re[n]=0;
      idft__XN.im[n]=0;
      for(k=0;k<N;k++)
       {
        //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
       angle=2*M_PI*k*n/N;
       idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
        idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
        }
     idft__XN.im[n]=idft__XN.im[n]/N;
      idft__XN.re[n]=idft__XN.re[n]/N;

     }
    */
     
    return idft__XN;
     }
     
     
     cplxvector  getanalyticsignal(double *xn, int num)
     {
      cplxvector  x_n,Hilb;
     int i;
       //x=x+jy,y=0
       //hilb(x)=a+jb
       //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
        //I=x(t), Q=j*Hilb(xre(t))

    Hilb=hilbert(xn,num);
     dcomp data1;
     double Ximg;
     for (i=0; i<num; i++)
      {
       data1=Hilb.at(i) ;
       Ximg=data1.imag();
       data1=(xn[i],Ximg );
       x_n.push_back(data1);
      //x_n.re[i]=xn[i];
      //x_n.im[i]=Hilb.im[i];
     }
     
     return x_n;
     }

     /*
     
    double Arg(double re, double im)
    {
     
    return   atan2 (im,re)  ;
    }
    */

    cplxvector  get_mod_fi(cplxvector x, int n, int mode   )
      {
      int i;
      cplxvector res;
      double fi;
      double module;
      for (i=0;i<n;i++)
      {
       module=abs(x.at(i));
       fi=arg(x.at(i));
      //res.re[i]=//pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
      //double treshh=1e9;
      // x.im[i]=x.im[i]*treshh;
      // x.im[i]=floor(x.im[i])/treshh;
      // x.re[i]=x.re[i]*treshh;
      // x.re[i]=floor(x.re[i])/treshh;
      //res.im[i]=Arg(x.re[i],x.im[i]);
       if (mode==1)  { fi=fi*180/M_PI;  }
        res.push_back((module,fi));
    std::cout<< "\n f="<<  " Hz;|Xf(jw)|=" << module << ";arg(Xf)="<< fi<<"  deg " ;
      }
       return  res;
     }

     
     
    cplxvector  GetSpectrum(double *xn, int num)
    {
     
    cplxvector  xn_iq, Xf_IQ, Xf;
    xn_iq=getanalyticsignal( xn,  num);
     
    // Xf_IQ=cdft(xn_iq  ,num );
     //Xf_IQ=ScaleDFT(Xf_IQ,num);
     // Xf =get_mod_fi(Xf_IQ, num,1   );
     
    return Xf;
    }


    /*


    Array getenvelope(double *x, int num)
     {
     int i;
      CData  y;
     Array env;
     y=hilbert(x,num);
      for (i=0;i<num;i++)
      {
      env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
      }
       return env;
     }
    */


    //FILE   *fp1,*fp2 ;


    int main()
    {

    double xin[NMAX];
    cplxvector  Xf;

    char key;
    int   n,i,N1;
    //int k;
    double fs,T, f;
    //double x;


    //printf ("\n For input  x(t) from file input '1' ");
    printf ("\n For manual  input  x(t) input   '0' ");
    std::cin>>key;
    //if (key=='1' ) { goto label_05;}
    if (key=='0' ) { goto label_06;}
    /*

    //label_05:

     printf("\n Open xinput.txt for reading x[n]   ");
    if ((fp1=fopen("xinput.txt","r"))==NULL)  printf ("\n File opening error ");
    fscanf(fp1,"%le",&fs);
    fscanf(fp1,"%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     fscanf (fp1,"%d",&k);
      fscanf (fp1,"%le",&xin[n]);
    }


    fclose(fp1);

    goto label_07;

    */

    label_06:

    printf("\n Input fs,Hz , ( %lf ) ",fs) ;
    std::cin>>fs; //scanf( "%le",&fs);
     printf("\n Input number of samples of xin[n]  ") ;
     std::cin>>N1; //scanf( "%d",&N1);

    for (n=0 ; n<N1;n++)
    {
     T=n*1/(N1*fs);
     std::cout<<"\n n="<<n<<"; t="<<T<<" s;   xin[n]= " ;
     std::cin>>xin[n];
     }

    // label_07:;

    std::cout<< "\n fs=" <<fs << " Hz " ;
    std::cout<< "\n Nsamples= "<<N1 ;

    for (n=0 ; n<N1;n++)
    {
     T=n*1/(N1*fs);
      std::cout<< "\n n="<< n << "; x[n]=" <<xin[n] << ",  t="<<T<<" s" ;
    }

    std::cout<<"\n Input 0 for continue " ;
    std::cin.get();
    std::cout<<"\n OK.   " ;


    /*****************************/

    //label1:;
    //Xf.re= |Xf|,  Xf.im -arg(Xf) (only in this program)
     std::cout<<"\n Obtaining spectrum of x[n]" ;
    Xf=GetSpectrum(xin, N1);

     std::cout<< "\n Spectrum of x[n] : " ;

     //if ((fp2=fopen("spectrum.txt","a"))==NULL)  printf ("\n File opening error");
     //fprintf(fp2,"\n   %le   %d \n",fs,N1);
     //fprintf(fp2, "\n n;\t f, Hz;\t |Xf(jw)| ;\t  arg(Xf), deg ");
    for (i=0; i<N1;i++)
    {
    f=fs*i/N1;
    //cout <<"\n f="<< f<<" Hz;|Xf(jw)|="<<Xf.re[i]<<";arg(Xf)="<<Xf.im[i]<<" deg ";
    // std::cout<< "\n f="<<f<<" Hz;|Xf(jw)|=" <<   real(Xf.at(i))<< ";arg(Xf)="<< imag(Xf.at(i))<<"  deg " ;
     //fprintf(fp2,"\n %d   %le  %le",i ,Xf.re[i],Xf.im[i]);

    }

    std::cout<< "\n  input   '0' for exit " ;
     
    std::cin.get();

    return 0;

    }

    Tuesday, February 6, 2018 6:25 PM
  • // RectAIMSpectr.cpp : 

    #include "stdafx.h"
    #include <stdio.h>
    #include <stdlib.h>
    #include <iostream>
    #include <math.h>

    using namespace std;
    #define M_PI 3.1415926535897932384626433832795


    double tau0, deltatau, A0;
    double f0, Fm,f1;
    double w0,T0,m,m1,Wm;

     


    int k,n;

    #define kmax 60 //change
    #define nmax 60 //change 

    double U0;
    double UFm[nmax], Ffm[nmax];
    double fi[kmax], Ufi[kmax];
    double fk[2*nmax+1][kmax+1], Ufk[2*nmax+1][kmax+1];
    double Q0;


    double Q;
    int kmax1;

    void InputData3()
    {


     //Q0=T0/tau0;
     //T0=1/f0
     //tau0=T0/Q0;

    std::cout<<"\nInput A0,V : (1) ";
    std::cin>>A0;
    std::cout<<"\nInput f1,Hz : (16384) ";
    std::cin>>f1;

    std::cout<<"\nInput Q  : (10) ";
    std::cin>>Q;
     std::cout<<"\nInput kmax  : (60) ";
    std::cin>>kmax1;
     
    return;
    }
    /*
     ________                    _______
    |             |                    |
    |             |___________|

    */
    void getFreqRectQ(void)
    {

    double U0=A0/Q;
    double timp=1/(f1*Q);
    double w1=2*M_PI*f1;
     
    double Ti=1/f1;
    std::cout <<"\n            f0="<<0<<"Hz ; U0=" <<U0<<" V  ";

    for(k=1;k<=kmax1; k++)
    {
    //Ufi~ A0*|sinc(k*pi/q)|  =A0*sin (k*pi/q) /(k*pi/q)
    //Ufi[k]=2*A0*sin(k*M_PI/Q)/(k*M_PI);
    //fi[k]=k*f1;

    //std::cout <<"\n k="<<k<<";  f[k]="<<fi[k]<<"Hz ; Um[k]=" <<Ufi[k]<<" V  ";

    double Um;
    //Um=fabs( ( 2*A0 /(k*M_PI))*sin(k*M_PI/Q));
     
    //std::cout <<"\n k="<<k<<";  f[k]="<<k*f1<<"Hz ; Um[k]=" <<Um<<" V  ";
     Um=(2*timp/Ti)*A0*(sin(k*2*M_PI*f1*timp/2)/(k*2*M_PI*f1*timp/2));
     Um=fabs(Um);
     std::cout <<"\n k="<<k<<";  f[k]="<<k*f1<<"Hz ; Um[k]=" <<Um<<" V  ";

    }

     std::cin.get() ;


    return;
    }

    /********************************************/
    /*

    S(t)=A0+sum(from k=1;to k=m;{ Ak*cos(k*wi*t)})=
    =(taui*Ui/Ti)+(2*taui/Ti*Ui)*sum(from k=1;to k=m;{
    (sin(k*wi*taui/2)/(k*wi*taui/2))*cos(kwit))

    uaim(t)=(taui*Ui/Ti)+(taui/Ti)*Ui*maim*cos(W*t)+
    +(2*taui/Ti)*Ui*sum(from k=1;to k=m;{
    (sin(k*wi*taui/2)/(k*wi*taui/2))*(1+maim*cos(W*t))*cos(k*wi*t) })

    Uk*cos(k*wi*t)+0.5*maim*Uk*cos(k*wi-W)*+0.5*maim*Uk*cos(k*wi+W)*t
    Uk=(2*taui/Ti)*Ui*(sin(k*wi*taui/2)/(k*wi*taui/2))

    f=0 Um[0]= taui*Ui/Ti
    f=W/2pi Um[W]=taui*miam*Ui/Ti

    f=k*wi-W Um=0.5*maim*Uk
    f=k*wi    Uk=Uk
    f=k*wi+W Um=0.5*maim*Uk
    ...
    Uenv_level[2*pi/taui]=0

    fsr=Whi

    */

    double maim;

    void InputData4()
    {


     //Q0=T0/tau0;
     //T0=1/f0
     //tau0=T0/Q0;

    std::cout<<"\nInput A0,V : (1) ";
    std::cin>>A0;
    std::cout<<"\nInput f1,Hz : (16384) ";
    std::cin>>f1;

    std::cout<<"\nInput Q  : (10) ";
    std::cin>>Q;
    std::cout<<"\nInput Fmod ,Hz : (1000) ";
    std::cin>>Fm;
    std::cout<<"\nInput maim,%  : (10) ";
    std::cin>>maim;
    maim/=100;
     std::cout<<"\nInput kmax  : (60) ";
    std::cin>>kmax1;
     
    return;
    }


    void getFreqRectAIM(void)
    {

     
    double timp=1/(f1*Q);
    double w1=2*M_PI*f1;
    double Ti=1/f1;

     
     
     
     double U0=timp*A0/Ti;
     std::cout <<"\n  ;  f[0]="<<0<<"Hz ; Um[k]=" << U0<<" V  ";
     double Fmod=Fm;
     double Ufm=timp*maim*A0/Ti;
     std::cout <<"\n  ;  f[Fm]="<<Fmod<<"Hz ; Um[k]=" << Ufm<<" V  ";
     double Uk,Ulsbk,Uusbk;
     double fk,flsbk,fusbk;
      std::cout <<"\n  ;  f[2*pi/taui  /( 2*pi)]="<< 1/timp<<"Hz ; Um[1/tau]=" << 0<<" V  ";
     for (int k=1;k<kmax1;k++){
     
     
     Uk=(2*timp/Ti)*A0*(sin(k*2*M_PI*f1*timp/2)/(k*2*M_PI*f1*timp/2));
     Uk=fabs(Uk);
     Ulsbk=0.5*maim*Uk;
     Uusbk=0.5*maim*Uk;
     fk=k*f1;
     flsbk=fk-Fm;
     fusbk=fk+Fm;
     
     std::cout <<"\n  ;  f[lsbk]="<<flsbk<<"Hz ; Um[flsbk]=" <<Ulsbk<<" V  ";
     std::cout <<"\n  ;  f[k]=   "<<fk<<" Hz   ; Um[fk]=" <<Uk<<" V  ";
     std::cout <<"\n  ;  f[usbk]="<<fusbk<<"Hz ; Um[fusbk]=" <<Uusbk<<" V  ";
     
     
     }

     std::cin.get() ;

    return;
    }

    int main( )
    {
     
     int key;
     std::cout<<"\n Input 1 for rect. signal\n Input 2 for AIM ";
     std::cin>>key;
     
     switch(key)
     {


     case 1:
             std::cout<<"\n  Rect. signal with Q   ";
            InputData3();
            getFreqRectQ();
           break;
     case 2:
         InputData4();
        getFreqRectAIM( );
       break;
     }

      std::cout<<"\n Input 0 for exit ";
      std::cin.get() ;
    return 0;
    }

    Tuesday, February 6, 2018 9:12 PM