-
Notifications
You must be signed in to change notification settings - Fork 3
/
tst_xfr.c
87 lines (73 loc) · 2.41 KB
/
tst_xfr.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
// tst_xfr.c - test transfer-function computaton
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "sigpro.h"
/**********************************************************/
void
freq_resp(float *H, int nf, double df, char *fn)
{
double f;
float *db, *ph, *gd;
int i;
FILE *fp;
db = (float *) calloc(nf, sizeof(float));
ph = (float *) calloc(nf, sizeof(float));
gd = (float *) calloc(nf, sizeof(float));
sp_cdb(H, db, nf); // decibels
sp_cph(H, ph, nf); // phase
sp_unwrap(ph, ph, nf); // unwrap phase
sp_cgd(H, gd, nf, df); // delay
fp = fopen(fn, "w");
if (!fp) {
printf("can't open %s\n", fn);
return;
}
fprintf(fp, "; %s - transfer-function test test\n\n", fn);
fprintf(fp, "; f dB ph gd\n");
for (i = 0; i < nf; i++) {
f = i * df;
fprintf(fp, "%8.3f %8.3f %8.3f %8.3f\n", f, db[i], ph[i], gd[i]);
}
fclose(fp);
free(db);
free(ph);
free(gd);
}
int
main()
{
double df, fc, sr;
float wn[2];
float *a, *b, *x, *y, *z;
int ft, no, np, nc, nf;
sr = 48000; // sampling rate
fc = 8000; // low-pass cut-off frequency
ft = 0; // filter type: 0=low_pass, 1=high_pass
no = 4; // filter order
np = 8192; // number of samples = buffer length
nc = no + 1; // number of filter coefficients
nf = np / 2 + 1; // number of frequencies
a = (float *) calloc(nc, sizeof(float));
b = (float *) calloc(nc, sizeof(float));
x = (float *) calloc(np, sizeof(float));
y = (float *) calloc(np, sizeof(float));
z = (float *) calloc(nf * 2, sizeof(float));
// compute white noise
sp_randflat(x, np); // noise stimulus in x
// Butterworth low-pass filter
wn[0] = (float)(2 * fc / sr); // cutoff_frequency / Nyquist_rate
sp_butter(b, a, no, wn, ft); // Butterworth filter coefficients
sp_filter(b, nc, a, nc, x, y, np); // filter response in y
// compute transfer function
sp_transfer(x, y, np, z); // transfer function in z
// write results
df = (sr / np) / 1000; // frequency resolution (kHz)
freq_resp(z, nf, df, "transfer.txt");
free(a);
free(b);
free(x);
free(y);
free(z);
return (0);
}