aboutsummaryrefslogtreecommitdiff
path: root/matlab/comms/script_verify_fir.m
blob: d53842f7af9aa1fb86757620951375ebb7a01fc4 (plain)
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
close all; clear *;

% Whether data shall be written to file or not
writeDataToFile = 0;

% Filter order
N = 7;
% Sampling frequency in <Hz>
fs = 50e6;
% Cutoff frequency in <Hz>
fc = 1e6;

% Design filter (floating-point)
h = fir1(N, fc/(fs/2));

% Plot frequency response
figure(1);
[H, f] = freqz(h, 1, 2^10, fs);
subplot(211);
plot(f/1e6, 20*log10(abs(H)), 'LineWidth', 2); grid on; hold on; box off;
xlabel('f in MHz \rightarrow'); ylabel('|H(f)| in dB \rightarrow');
subplot(212);
plot(f/1e6, unwrap(angle(H)), 'LineWidth', 2); grid on; hold on; box off;
xlabel('f in MHz \rightarrow'); ylabel('arg H(f) in rad \rightarrow');

% Plot impulse response
figure(2);
stem(0:N, h, 'LineWidth', 2); hold on;
xlabel('Coefficient index \rightarrow'); ylabel('Amplitude \rightarrow');

% Wordlength of coefficients
w_c = 12;

% Convert floating-point coefficients to fixed-point
hqi = round(h*2^(w_c-1));
hq = hqi/2^(w_c-1);

% Save integer coefficients to file for later use in VHDL model
if writeDataToFile == 1
  fileID = fopen('hqi.txt','w');
  fprintf(fileID, ['CONSTANT COEFFS : COEFFS_TYPE := (']);
  for k = 1:N
    fprintf(fileID, '%d,', hqi(k));
  end
  fprintf(fileID, '%d);\n', hqi(N+1));
  fclose(fileID);
end

% Plot quantized frequency response
figure(1);
[Hq, f] = freqz(hq, 1, 2^10, fs);
subplot(211);
plot(f/1e6, 20*log10(abs(Hq)), 'LineWidth', 2);
title(['Transfer function, N = ' num2str(N), ', fc = ' num2str(fc/1e6) ' MHz']);
subplot(212);
plot(f/1e6, unwrap(angle(Hq)), 'LineWidth', 2);

% Plot quantized impulse response
figure(2);
stem(0:N, hq, 'x', 'LineWidth', 2); hold off;
legend('Float', ['Fixed 1.' num2str(w_c-1) 's, normalized']);
title('Impulse response');

% Create input sample stream
nofSamples = 2^12;
sigma = 0.25;
x = sigma*randn(1, nofSamples);

% Find outliers and saturate to [-1 1-LSB]
w_in = 14;
idx = (x >= 1); x(idx) = 1-2^-(w_in-1);
idx = (x < -1); x(idx) = -1;

% Display histogram
figure(3);
histogram(x, 100);
xlabel('x \rightarrow'); ylabel('Occurences \rightarrow');
title(['Distribution of input samples, \sigma = ' num2str(sigma)]);

% Quantize input sample stream to w_in bit
xqi = round(x*2^(w_in-1));
xq = xqi/2^(w_in-1);

% Save input sample stream to file
if writeDataToFile == 1
  save_variable(xqi, '%d', '../../sim/fir/stimuli/fir_stimuli.dat');
end

% Filter quantized input stream with quantized impulse response
yq = filter(hq, 1, xq);

% Plot quantized input and output sample streams
figure(4);
t = (0:nofSamples-1)/fs;
plot(t*1e6, xq); hold on; box off;
plot(t*1e6, yq);
xlabel('t in microseconds \rightarrow'); ylabel('Amplitude \rightarrow');

% Estimate transfer function given specific outcome of xq
[Hxyq, f] = tfestimate(xq, yq, [], [], 2^10, fs);

% Plot estimated frequency response
figure(1);
subplot(211);
plot(f/1e6, 20*log10(abs(Hxyq)), 'LineWidth', 2);
subplot(212);
plot(f/1e6, unwrap(angle(Hxyq)), 'LineWidth', 2);

% Truncation of output result to w_out bit
w_out = 14;

% Filter quantized input with quantized impulse response
yqi = filter(hqi, 1, xqi);

% Determine dynamic wordlength of current result
w_dyn_cur = ceil(log2(max(abs(yqi)))) + 1;
disp(['Dynamic wordlength, current result: ' num2str(w_dyn_cur) ' bits']);

% Determine worst-case dynamic wordlength
w_dyn_wc = w_in + floor(log2(sum(abs(hqi)))) + 1;
disp(['Dynamic wordlength, worst-case: ' num2str(w_dyn_wc) ' bits']);
yqit = floor(yqi/2^(w_dyn_wc-w_out));

% Determine dynamic wordlength of current output
w_dyn_curs = ceil(log2(max(abs(yqit)))) + 1;
disp(['Dynamic wordlength, current output: ' num2str(w_dyn_curs) ' bits']);

% Plot truncated output sample stream
figure(4);
plot(t, yqit / 2^(w_out-1));
legend('x(t)', 'y(t)', 'y(t), truncated');
title('Quantized input and output samples');

% Estimate transfer function given specific outcome of xqint
[Hxyqints, f] = tfestimate(xqi, yqit, [], [], 2^10, fs);

% Plot the estimated frequency response
figure(1);
subplot(211);
plot(f/1e6, 20*log10(abs(Hxyqints)*2^(w_in-w_out)), 'LineWidth', 2);
legend('Float', ['Fixed 1.' num2str(w_c-1) 's'], 'Estimate', 'Estimate, truncated', 'Location','southwest');
subplot(212);
plot(f/1e6, unwrap(angle(Hxyqints)), 'LineWidth', 2);

% Save output sample stream to file
if writeDataToFile == 1
  save_variable(yqit, '%d', '../../sim/fir/log/fir_result_ref.dat');
end