2 files changed, 164 insertions, 0 deletions
diff --git a/matlab/comms/save_variable.m b/matlab/comms/save_variable.m
new file mode 100755
index 0000000..6bc45e6
--- /dev/null
+++ b/matlab/comms/save_variable.m
@@ -0,0 +1,16 @@
+function save_variable (var, format, filename)
+
+% function save_variable (var, format, filename)
+%
+% Saves a variable var to filename using format,
+% e.g., save_variable (x, '%d', 'input.dat');
+
+[fid, message] = fopen(filename, 'w');
+if fid == -1
+	disp('File error. Message returned was:')
+	disp(message)
+	return
+end
+fprintf(fid, [format,'\n'], var);
+fclose(fid);
+return;
diff --git a/matlab/comms/script_verify_fir.m b/matlab/comms/script_verify_fir.m
new file mode 100755
index 0000000..d53842f
--- /dev/null
+++ b/matlab/comms/script_verify_fir.m
@@ -0,0 +1,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
+\ No newline at end of file