From d7b4d9beb249c92f61fb28a9fece02aeb084e24f Mon Sep 17 00:00:00 2001
From: Ksyer <>
Date: Thu, 9 May 2024 16:48:29 +0800
Subject: [PATCH] Add experiment

---
 1 - experiment/Chirp/Test.m                   | 20 ++++++++
 1 - experiment/Chirp/compress.m               | 17 +++++++
 1 - experiment/Chirp/configure_parameters.m   | 21 +++++++++
 1 - experiment/Chirp/frequency_domain.m       | 39 +++++++++++++++
 1 - experiment/Chirp/time_domain.m            | 24 ++++++++++
 .../Circulant Matrix/circulant_matrix.m       | 47 +++++++++++++++++++
 .../Circulant Matrix/is_vector_equal.m        |  5 ++
 7 files changed, 173 insertions(+)
 create mode 100644 1 - experiment/Chirp/Test.m
 create mode 100644 1 - experiment/Chirp/compress.m
 create mode 100644 1 - experiment/Chirp/configure_parameters.m
 create mode 100644 1 - experiment/Chirp/frequency_domain.m
 create mode 100644 1 - experiment/Chirp/time_domain.m
 create mode 100644 1 - experiment/Circulant Matrix/circulant_matrix.m
 create mode 100644 1 - experiment/Circulant Matrix/is_vector_equal.m

diff --git a/1 - experiment/Chirp/Test.m b/1 - experiment/Chirp/Test.m
new file mode 100644
index 0000000..2ba5e52
--- /dev/null
+++ b/1 - experiment/Chirp/Test.m	
@@ -0,0 +1,20 @@
+N = 1e3 + 1;
+f_c = 1e6;
+t = linspace(0, pi / f_c, N);
+f_s = 1 / N;
+freqs = (0:N-1) / f_s;
+x = exp(-1j * 2 * pi * f_c * t);
+
+subplot(211);
+plot(t, imag(x));
+title("x(t) = sin(t)");
+xlabel("t");
+ylabel("x(t)");
+
+x_fft = abs(fft(imag(x)));
+subplot(212);
+plot(t, x_fft);
+title("X(f) = fft(x(t))");
+xlabel("f");
+ylabel("X(f)");
+
diff --git a/1 - experiment/Chirp/compress.m b/1 - experiment/Chirp/compress.m
new file mode 100644
index 0000000..b9a4355
--- /dev/null
+++ b/1 - experiment/Chirp/compress.m	
@@ -0,0 +1,17 @@
+clc; clear;
+
+configure_parameters;
+global s_T s_R
+
+compressed_signal = conv(s_R, fliplr(s_T), 'same');
+
+figure(1);
+subplot(311);
+plot(t, s_T);
+title("Chirp signal");
+subplot(312);
+plot(t, s_R);
+title("Received Chirp signal");
+subplot(313);
+plot(t, compressed_signal);
+title("Matcher filter");
diff --git a/1 - experiment/Chirp/configure_parameters.m b/1 - experiment/Chirp/configure_parameters.m
new file mode 100644
index 0000000..cbee333
--- /dev/null
+++ b/1 - experiment/Chirp/configure_parameters.m	
@@ -0,0 +1,21 @@
+global f_s T_r simu_time t start_freq end_freq c len r_0 freqs s_T s_R
+
+f_s = 1e4;
+T_r = 1 / f_s;
+simu_time = 1;
+t = linspace(0, simu_time, f_s);
+start_freq = 0;
+end_freq = 1e7;
+c = 3e8;
+r_0 = 8e6;
+
+s_T = chirp(t, start_freq, simu_time, end_freq, 'linear') * 10;
+len = length(s_T);
+freqs = linspace(-f_s / 2, f_s / 2, len);
+
+s_R = zeros(size(s_T));
+delay = 2 * r_0 / c;
+delay_samples = round(delay * len) + 1;
+s_R(delay_samples:len) = 0.8 * s_T(1:len-delay_samples + 1);
+s_R = awgn(s_R, 1);
+
diff --git a/1 - experiment/Chirp/frequency_domain.m b/1 - experiment/Chirp/frequency_domain.m
new file mode 100644
index 0000000..7252103
--- /dev/null
+++ b/1 - experiment/Chirp/frequency_domain.m	
@@ -0,0 +1,39 @@
+clc; clear;
+
+configure_parameters;
+global s_T s_R c
+
+% 频域上的匹配滤波
+s_T_fft = fft(s_T);
+s_R_fft = fft(s_R);
+s_R_fft_conj = conj(s_R_fft);
+
+matched_filter = ifft(s_T_fft .* s_R_fft_conj);
+[~, peak] = max(matched_filter);
+t_delay = 1 - peak / length(matched_filter);
+d = (c * t_delay) / 2;
+
+fprintf("r = %f, predict d = %f\n", r_0, d);
+fprintf("eps = %f\n", abs(r_0 - d));
+
+figure(1);
+subplot(311);
+plot(freqs, abs(s_T_fft));
+title("Chirp signal");
+subplot(312);
+plot(freqs, abs(s_R_fft));
+title("Received Chirp signal");
+subplot(313);
+plot(freqs, abs(s_T_fft .* s_R_fft_conj));
+title("Matcher filter");
+
+figure(2);
+subplot(311);
+plot(t, s_T);
+title("Chirp signal");
+subplot(312);
+plot(t, s_R);
+title("Received Chirp signal");
+subplot(313);
+plot(t, matched_filter);
+title("After matched filter");
diff --git a/1 - experiment/Chirp/time_domain.m b/1 - experiment/Chirp/time_domain.m
new file mode 100644
index 0000000..08053e5
--- /dev/null
+++ b/1 - experiment/Chirp/time_domain.m	
@@ -0,0 +1,24 @@
+clc; clear;
+
+configure_parameters;
+global s_T s_R c
+
+% 时域上的匹配滤波
+matched_filter = conv(s_R, fliplr(s_T), 'same');
+[~, peak] = max(matched_filter);
+t_delay = peak / length(matched_filter) - 0.5;
+d = (c * t_delay) / 2;
+
+fprintf("r = %f, predict d = %f\n", r_0, d);
+fprintf("eps = %f\n", abs(r_0 - d));
+
+figure(1);
+subplot(311);
+plot(t, s_T);
+title("Chirp signal");
+subplot(312);
+plot(t, s_R);
+title("Received Chirp signal");
+subplot(313);
+plot(t, matched_filter);
+title("Matcher filter");
diff --git a/1 - experiment/Circulant Matrix/circulant_matrix.m b/1 - experiment/Circulant Matrix/circulant_matrix.m
new file mode 100644
index 0000000..f7837de
--- /dev/null
+++ b/1 - experiment/Circulant Matrix/circulant_matrix.m	
@@ -0,0 +1,47 @@
+%%
+clc;
+clear;
+
+N = 100;
+eps = 1e-10;
+F_n = dftmtx(N);
+
+for i = 1: N
+    for j = 1: N
+        fprintf(string(F_n(i, j)) + ", ");
+    end
+    fprintf("\b\b\n");
+end
+
+fprintf("\n");
+
+% 检查傅里叶矩阵的正确性
+assert(abs(F_n(2, 2) ^ 2 - F_n(2, 3)) < eps);
+assert(abs(F_n(2, 2) ^ N - 1) < eps);
+
+%% 检查循环卷积
+x = randn(N, 1);
+z = randn(N, 1);
+p = cconv(x, z, N);
+
+% 检查循环卷积公式的正确性
+for i = 1: N
+    tmp = 0;
+    for j = 1: N
+        tmp = tmp + x(j) * z(mod(i-j+N, N)+1);
+    end
+    assert(abs(tmp - p(i)) < eps);
+end
+
+% 检查 x * z = A(x)z 的正确性
+A_x = toeplitz([x(1) fliplr(x(2:end)')], x);
+assert(is_vector_equal(p, A_x' * z, eps));
+
+% 检查卷积定理的正确性
+x_hat = fft(x);
+z_hat = fft(z);
+xz_hat = fft(p);
+assert(is_vector_equal(x_hat .* z_hat, xz_hat, eps));
+
+
+
diff --git a/1 - experiment/Circulant Matrix/is_vector_equal.m b/1 - experiment/Circulant Matrix/is_vector_equal.m
new file mode 100644
index 0000000..cf3cac1
--- /dev/null
+++ b/1 - experiment/Circulant Matrix/is_vector_equal.m	
@@ -0,0 +1,5 @@
+function flg = is_vector_equal(A, B, eps)   
+    diff = abs(A - B);
+    max_diff = max(diff);
+    flg = max_diff <= eps;
+end
\ No newline at end of file