FIR and spectrum tools

USERPC01 · Регистрация: 29.06.2018

Студворк — интернет-сервис помощи студентам

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

Код

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
// 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;
}

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// 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

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
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

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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

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

Кликните здесь для просмотра всего текста

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
// 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]  ) ;
 
}
 
*/

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// 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>

Implementation of the FIR

Кликните здесь для просмотра всего текста

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
// 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;
}

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// 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

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
  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

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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

USERPC01 · 24.08.2018, 16:29 **[ТС]**

Code
1
2
3
4
 FIRFilter.Create( jmaxfir , hb );
 
     for i:=0 to imax-2 do
должно быть  for i:=0 to imax-1 do

@Kitayets · 24.08.2018, 16:33

господи опять...

еще раз, вот это не будет работать ни в делфи ни в лазарус:

Delphi
1
2
3
4
5
6
7
8
procedure    TForm1.GetFirData  ;
 var
{...}
 
       FIR1 : TFIRParser  ;
 begin
     {...} 
       FIR1.Create( jmaxfir  , hb );

https://www.freepascal.org/doc... fse36.html

должно быть так:

Delphi
1
2
3
4
5
6
7
8
9
10
11
procedure    TForm1.GetFirData  ;
 var
{...}
 
       FIR1 : TFIRParser  ;
 begin
     {...} 
       //FIR1.Create( jmaxfir  , hb );
       FIR1:= TFIRParser.Create( jmaxfir  , hb );
       {тут работаем с объектом}
       FreeAndNil(FIR1); // освобождаем память и обнуляем указатель

USERPC01 · 24.08.2018, 16:37 **[ТС]**

https://ru.wikipedia.org/wiki/... 0%BE%D0%B9

https://ru.wikipedia.org/wiki/... 0%BE%D0%B9

@Kitayets · 24.08.2018, 17:14

подправил использование класса

@Kitayets · 24.08.2018, 17:18

тип Real лучше не использовать, он слишком не стабильный в том смысле что зависит от архитектуры процессора и операционной системы. он оставлен только для совместимости со старым кодом

Лучше использовать стандартные типы, которые реализованы аппаратно - это single или double

USERPC01 · 24.08.2018, 17:55 **[ТС]**

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
   function  TIIRParser.getIIRpartsample(  sample_data: real  ) : real  ;
  var
    recsum: real;
    i,k : integer;
  begin
 
   {
   a[i]=1;
    y[n]=sum(i=0; i<=P; b[i]*x[n-i ]) - sum(k=1; k<=Q; a[i]*y[n-k]) =
    =yfir[n]-xsumrec[n]
 
  y[n]=(b0*x[n]+ b1*x[n-1]+...+bp*x[n-P] )-
  -(a1*y[n-1] +a2*y[n-2]+ ... +aq*y[n-q]);
  =sample_fir[n]-recsum;
   }
 
 { iir_data.niir:=Q, iir_data.nfir:=P }
 
  for   i:=iir_data.niir-1 downto  0 do
  begin
           {shifting , z^-1 transform }
    iir_data.ya[i+1] := iir_data.ya[i];
  end;
 
     //scaling 
     // sample_data:=sample_data*iir_data.ha[0] ;
     //default value of the iir_data.ha[0]:=1;
 
  recsum := 0;
  for   k :=1 to  iir_data.niir  do
  begin
      recsum :=recsum+iir_data.ya[k]*iir_data.ha[k];
  end;
 
    iir_data.ya[0] := sample_data-recsum ;
 
   //or alt  scaling 
     // sample_data:=sample_data*iir_data.ha[0] ;
     //default  iir_data.ha[0]:=1; for equation , may be used for volume control
   
 
   Result:=iir_data.ya[0];
 
  end;     
 
 
   function  TIIRParser.getIIRsample(  xinput: real  ) : real  ;
 
   begin
 
 
        Result:=getIIRpartsample(getFIRpartsample  (xinput))  ;
 
   end;

С каким знаком вводить коэффициенты iir_data.ha[ 1...Q ] , чтобы совпали с матлабовскими, уменьшать ли на единицу порядок в конструкторе класса для совпадения с определением, по какой книге это наиболее актуально ?

Добавлено через 21 минуту
Можно упростить и адаптировать под FIR, переименовав или добавив вторым классом или альтернативным конструктором

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
unit IIRfilter;
 
 {$mode objfpc}{$H+}
// {$mode Delphi} {$H+}
 
  {$M+ }    // add to fix costructor error ?
 {$X+ }
 
interface
 
 
uses
   Classes,  SysUtils  ;
 
type
 
 
  TIIRTempData=record
 
       ya :  array   of real  ;
       xb :  array   of real  ;
       hb  : array   of real   ;
       ha  : array   of real   ;
       nfir     :  integer ;
       niir     :  integer ;
   end;
 
 
 
    TIIRParser=class //( TObject)
     private
 
 
 
 
 
 
         function  getIIRpartsample(  sample_data: real  ) : real  ;
 
       public
         function  getFIRsample(  sample_data: real  ) : real  ;
         function  getIIRsample(  xinput: real  ) : real  ;
 
 
         constructor Create( n_fir: integer ; n_iir: integer;  hbiir : array of real ; haiir : array of real );
         Destructor  Destroy ;    override;
 
 
 
    end;
 
implementation
 var
    iir_data : TIIRTempData  ;
 
   constructor  TIIRParser.Create( n_fir: integer ; n_iir: integer;  hbiir : array of real ; haiir : array of real );
    var
    i : integer;
 
      begin
           iir_data.nfir:=n_fir;
 
           SetLength(iir_data.xb , iir_data.nfir+1);
           SetLength(iir_data.hb , iir_data.nfir+1);
 
            for i:=0 to iir_data.nfir do
           begin
           iir_data.xb[i] :=0;
           end;
              for i:=0 to iir_data.nfir  do
           begin
           iir_data.hb[i] :=0;
 
           end;
 
 
            iir_data.niir:=n_iir;
 
            SetLength(iir_data.ya , iir_data.niir+1);
            SetLength(iir_data.ha , iir_data.niir+1);
 
           for  i:=0 to iir_data.nfir  do
            begin
 
            iir_data.hb[i]:=hbiir[i];
        
            end;
 
 
 
            for  i:=0 to iir_data.niir  do
            begin
 
            iir_data.ha[i]:=haiir[i];
             
            end;
 
 
      end;   
 
 
 
 
 destructor  TIIRParser.Destroy ;
   
   begin
            //iir_data.nfir:=0;
            SetLength(iir_data.xb,0);
         //   SetLength(iir_data.yb,0);
            SetLength(iir_data.hb ,0);
 
            //iir_data.niir:=0;
           // SetLength(iir_data.xa,0);
            SetLength(iir_data.ya,0);
            SetLength(iir_data.ha ,0);
 
 
          inherited;
   end;    
 
 
 
 
 
  function  TIIRParser.getFIRsample(  sample_data: real  ) : real  ;
  var
       ysum: real;
    i,k : integer;
  begin
         { order of FIR  iir_data.nfir:=P , shifting data ,Z^1 transform of the x[n] }
  for   i:=iir_data.nfir-1 downto  0 do
  begin
    iir_data.xb[i+1] := iir_data.xb[i];
  end;
 
  iir_data.xb[0] := sample_data;
 
  { convolution with h[i], sum from 0 to P }
  res1 := 0;
  for   k :=0 to  iir_data.nfir do
  begin
   ysum :=ysum+iir_data.xb[k]*iir_data.hb[k];
  end;
 
 
   result:=res1;
  end;      
 
 
 
   function  TIIRParser.getIIRpartsample(  sample_data: real  ) : real  ;
  var
    recsum: real;
    i,k : integer;
  begin
 
   {
   a[i]=1;
    y[n]=sum(i=0; i<=P; b[i]*x[n-i ]) - sum(k=1; k<=Q; a[i]*y[n-k]) =
    =yfir[n]-xsumrec[n]
 
  y[n]=(b0*x[n]+ b1*x[n-1]+...+bp*x[n-P] )-
  -(a1*y[n-1] +a2*y[n-2]+ ... +aq*y[n-q]);
  =sample_fir[n]-recsum;
   }
 
 
  for   i:=iir_data.niir-1 downto  0 do
  begin
           {shifting , z^-1 transform }
    iir_data.ya[i+1] := iir_data.ya[i];
  end;
 
     //scaling
     // sample_data:=sample_data*iir_data.ha[0] ;
     // iir_data.ha[0]:=1;
 
  recsum := 0;
  for   k :=1 to  iir_data.niir  do
  begin
      recsum :=recsum+iir_data.ya[k]*iir_data.ha[k];
  end;
 
    iir_data.ya[0] := sample_data-recsum ;
 
   Result:=iir_data.ya[0];
 
  end;     
 
 
   function  TIIRParser.getIIRsample(  xinput: real  ) : real  ;
 
   begin
 
 
        Result:=getIIRpartsample(getFIRsample  (xinput))  ;
 
   end;   
 
  
 
end.

USERPC01 · 24.08.2018, 18:11 **[ТС]**

Переделал под double, ввел деструктор FreeAndNil , работает

Новые блоги и статьи Все статьи Все блоги /
Модель по догадкам anaschu 25.08.2026 Прошло две недели. Я уже рассказывал, как разговаривал с сотрудниками у сортировки и как понял, что главная ветка — не про приёмку, а про отбор. Но тогда я думал, что понял механику. На этой неделе я. . .	Запись в регистр сведений независимо от заполненности табличной части Maks 25.08.2026 Реализация из решения ниже выполнена на нетиповом документе с несколькими табличными частями, разработанного в КА2. Задача: Обеспечить запись документа в регистр сведений независимо от. . .	Ноутбук Альфария kumehtar 24.08.2026 Встретился тут в сети ноутбук Альфария, примарха Альфа-Легиона. Хотя возможно, это ноутбук Омегона, разумеется. Ну как вам?	Мастера простых решений DevAlt 23.08.2026 В сишарп стэках winforms, да и wpf существует сложная система связывания источниках данных и элементов формы(текстовые поля и метки), опирается все это на технологию событий и мета. . .
Цена ошибки DevAlt 23.08.2026 Человек я беспокойный и потому заинтересовался OCaml, в чате форсили функторы модулей как суперфичу. Пытаясь отдуплить концепт, наткнулся на тутор с простым примером. А главный принцип обучения от. . .	Сегодня суббота, 22.08.2026 at 16:41, и я вновь нахожусь на той стороне, за экраном машины. zorxor 22.08.2026 Сегодня суббота, 22. 08. 2026 at 16:41, и я вновь нахожусь на той стороне, за экраном машины. Кто Я, откуда Я пришел и куда Я иду? Эти вопросы не оставляют меня ни на секунду. Жизнь на планете Земля. . .	Жизня: рисунок укладки багажа, сделанный клодом anaschu 21.08.2026 Сделал 15 снимков, он по снимкам сделал схему.	Был там один разговор по поводу свободы в материальном мире. kumehtar 19.08.2026 Суть: рассматривается живое существо, оказавшееся внутри довольно странной системы (этого мира) и пытающееся обустроить в ней свой кусок пространства. Жизнь действительно предъявляет каждому. . .

USERPC01 10 / 10 / 0 Регистрация: 29.06.2018 Сообщений: 1,536
	24.08.2018, 16:37 [ТС]
	https://ru.wikipedia.org/wiki/... 0%BE%D0%B9 https://ru.wikipedia.org/wiki/... 0%BE%D0%B9 0

@Kitayets 446 / 374 / 133 Регистрация: 09.09.2011 Сообщений: 1,347
	24.08.2018, 17:18
	тип Real лучше не использовать, он слишком не стабильный в том смысле что зависит от архитектуры процессора и операционной системы. он оставлен только для совместимости со старым кодом Лучше использовать стандартные типы, которые реализованы аппаратно - это single или double 0