FIR and spectrum tools

// 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

// 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>

// 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