信原matlab实验 (PKUanonym#51)

大三上/信号处理原理/hw/matlab/2019/ex1/ex1.log

@@ -0,0 +1,20 @@
+Single button test
+|length |fft-result |fft-time |grt-result |grt-time |time ratio |
+| ----  |   ----    |   ----  |   ----    |   ----  |    ----   |
+|25695  |     0     |0.002852 |     0     |0.001258 | 2.267032  |
+|14988  |     1     |0.000599 |     1     |0.000753 | 0.796466  |
+|19985  |     2     |0.004619 |     2     |0.001288 | 3.586335  |
+|19271  |     3     |0.001731 |     3     |0.000950 | 1.822278  |
+|19985  |     4     |0.001775 |     4     |0.001007 | 1.762288  |
+|23018  |     5     |0.002294 |     5     |0.001170 | 1.959672  |
+|19271  |     6     |0.001924 |     6     |0.000983 | 1.956478  |
+|23018  |     7     |0.000978 |     7     |0.001117 | 0.875862  |
+|18736  |     8     |0.000715 |     8     |0.001087 | 0.657708  |
+|22483  |     9     |0.002365 |     9     |0.001126 | 2.100906  |
+|23019  |     *     |0.002585 |     *     |0.001172 | 2.205220  |
+|24625  |     #     |0.002101 |     #     |0.001710 | 1.228426  |
+telephone number
+FFT time_cost: 0.034963 s
+ans: ---8-22-1177-5-6699-0-99-99
+Goertzel time_cost: 0.021653 s
+ans: 1--8-22-1177-5-6699-0099-99

大三上/信号处理原理/hw/matlab/2019/ex1/main.m

@@ -0,0 +1,125 @@
+% 主函数
+function main
+    global keys; % 按键表
+    global freqs; % 频率表
+    keys = ['1', '2', '3', 'A';
+            '4', '5', '6', 'B';
+            '7', '8', '9', 'C';
+            '*', '0', '#', 'D';
+            '-', '-', '-', '-'];
+    freqs = [697, 770, 852, 941, 1209, 1336, 1477, 1633];
+
+    % 输出到文件 ex1.log
+    fout = fopen('./ex1.log', 'w');
+
+    % 单按键测试
+    fprintf(fout, "Single button test\n");
+    fprintf(fout, "|length |fft-result |fft-time |grt-result |grt-time |time ratio |\n");
+    fprintf(fout, "| ----  |   ----    |   ----  |   ----    |   ----  |    ----   |\n");
+    for i = 0:11
+    	[audio, fs] = audioread(['./data/',num2str(i),'.wav']);
+        audio = sum(audio, 2) / 2; % 合并左右声道
+        % audio = audio(1:4000);
+    	fprintf(fout, "|%5d  ", length(audio));
+    	tic
+        ans = fft_test(audio, fs);
+        fft_time = toc;
+        fprintf(fout, "|     %c     |%f ", ans, fft_time);
+        tic
+        ans = goertzel(audio, fs);
+        grt_time = toc;
+        fprintf(fout, "|     %c     |%f ", ans, grt_time);
+        fprintf(fout, "| %f  |", fft_time/grt_time);
+        fprintf(fout, "\n");
+    end
+
+    % 一组电话号码测试
+    [audio, fs] = audioread('./data/18217569099.wav');
+    patch_size = fs/3;
+    audio = sum(audio, 2) / 2;
+    total_len = size(audio, 1);
+
+    fprintf(fout, "telephone number\n");
+    fft_result = [];
+    tic
+    for i = 1:patch_size:total_len-patch_size
+    	fft_result = [fft_result, fft_test(audio(i:i+patch_size), fs)];
+    end;
+    fprintf(fout, "FFT time_cost: %f s\n", toc);
+    fprintf(fout, "ans: %s\n", fft_result);
+
+    grt_result = [];
+    tic
+    for i = 1:patch_size:total_len-patch_size
+    	grt_result = [grt_result, goertzel(audio(i:i+patch_size), fs)];
+    end;
+    fprintf(fout, "Goertzel time_cost: %f s\n", toc);
+    fprintf(fout, "ans: %s\n", grt_result);
+
+    fclose(fout);
+end
+
+% FFT Algorithm
+function ans = fft_test(audio, fs)
+    global keys;
+    global freqs;
+
+    y = abs(fft(audio));
+    len = length(y);
+    % 数组下标是数字频率，需要转换成模拟频率
+    l1 = round(650  / fs * len);
+    h1 = round(1000 / fs * len);
+    l2 = round(1150 / fs * len);
+    h2 = round(1700 / fs * len);
+    [px, col] = max(y(l2:h2));
+    [py, row] = max(y(l1:h1));
+    col = col + l2 - 1;
+    row = row + l1 - 1;
+    fx = col / len * fs;
+    fy = row / len * fs;
+
+    ans = '-';
+    if (px > 10 || py > 10) % 过滤噪声
+        [t, row] = min(abs(freqs(1:4) - fy));
+        [t, col] = min(abs(freqs(5:8) - fx));
+        ans = keys(row, col);
+    end;
+end
+
+% Goertzel Algorithm
+function ans = goertzel(audio, fs)
+    global keys;
+    global freqs;
+
+    t_freqs = freqs;
+    N = length(audio);
+    x = audio;
+    v = zeros(1,N+2);
+    y = zeros(1, 8);
+    t_freqs = round(t_freqs * N / fs);
+
+    % 差分方程
+    % v_k[n] = 2cos(2k\pi/N)v_k[n-1] - v_k[n-2] + x[n]
+    % y_k[n]^2 = v_k[n] - W_N^k*v_k[n-1]
+    for k = 1 : length(t_freqs)
+        factor = 2 * cos(2 * pi * t_freqs(k) / N);
+        for i = 1 : N
+            v(i+2) = factor*v(i+1) - v(i) + x(i);          % v_k[n]
+        end
+        y(k) = v(N+2)^2 + v(N+1)^2 - factor*v(N+2)*v(N+1); % |y_k[N]|^2
+    end
+    y = sqrt(abs(y));
+
+    % 过滤噪声
+    [t, col] = max(y(5:8));
+    [t, row] = max(y(1:4));
+    px = y(col+4);
+    py = y(row);
+    trhld = max(px, py)*0.3;
+    if (sum((y-trhld) > 0) > 2)
+        row = 5;
+    end
+
+    ans = keys(row, col);
+end
+

大三上/信号处理原理/hw/matlab/2019/ex2/main.m

@@ -0,0 +1,98 @@
+function main
+	points = 400;
+	step = 100;
+	M = 1000; % length y
+	time_origin = zeros(1, points);
+	time_circle = zeros(1, points);
+	time_save   = zeros(1, points);
+	time_add    = zeros(1, points);
+
+	for i = 1 : points
+		N = M + (i-1) * step
+		fprintf("testing N = %d, M = %d\n", N, M);
+		x = rand([1, N]);
+		y = rand([1, M]);
+		tic;
+		ans1 = conv_origin(x, y);
+		time_origin(i) = toc;
+		tic;
+		ans2 = conv_circle(x, y);
+		time_circle(i) = toc;
+		tic;
+		ans3 = overlap_save(x, y);
+		time_save(i) = toc;
+		tic;
+		ans4 = overlap_add(x, y);
+		time_add(i) = toc;
+	end
+	x = M:step:M + step*(points-1);
+
+	plot(x, time_origin, '-', x, time_circle, '-', x, time_save, '-', x, time_add, '-');
+	legend({'origin conv', 'circle conv', 'overlap-save', 'overlap-add'}, 'Location', 'northwest');
+	xlabel('M');
+	ylabel('time_cost/s');
+end
+
+function ans = conv_origin(x, y)
+	N = length(x);
+	M = length(y);
+	L = M + N - 1;
+	ans = zeros(1, L);
+
+	for n = 1 : L
+		m0 = max(1, n + 1 - N);
+		m1 = min(M, n);
+		for m = m0 : m1
+			ans(n) = ans(n) + x(n - m + 1) * y(m);
+		end
+	end
+end
+
+function ans = conv_circle(x, y)
+	N = length(x);
+	M = length(y);
+	x = [x, zeros(1, M-1)]; % zero padding
+	y = [y, zeros(1, N-1)]; % zero padding
+	ans = ifft(fft(x) .* fft(y));
+end
+
+function y = overlap_add(x, h)
+    N = length(x);
+    M = length(h);
+    L = M + 1; % block size
+    len = L + M - 1;
+    fft_h = fft(h, len);
+    num = ceil(N / L); % split x into $num blocks
+    y = zeros(1, (num + 1) * L);
+    x = [x, zeros(1, (num + 1) * L - N)]; % zero padding to (num + 1) * L
+    ol = zeros(1, M - 1); % overlap-add
+    for i = 1:L:num * L + 1
+        y_k = ifft(fft(x(i:i + L - 1), len) .* fft_h); % calc circle conv
+        y_k(1:M - 1) = y_k(1:M - 1) + ol(1:M - 1); % overlap-add
+        ol(1:M - 1) = y_k(L + 1:len); % refresh
+        y(i:i + L - 1) = y_k(1:L); % save y
+    end
+    y = y(1:N + M - 1);
+end
+
+% conv by overlap-save
+function y = overlap_save(x, h)
+    % assume h is much shorter than x
+    N = length(x);
+    M = length(h);
+    L = M + 1; % block size
+    len = L + M - 1;
+    fft_h = fft(h, len);
+    num = ceil(N / L); % split x into $num blocks
+    y = zeros(1, (num + 1) * L);
+    x = [x, zeros(1, (num + 1) * L - N)]; % zero padding
+    ol = zeros(1, M - 1); % overlap-save
+    for i = 1:L:num * L + 1
+        subx = [ol, x(i:i + L - 1)];
+        ol = subx(L + 1 : L + M - 1); % refresh overlap-save
+        y_k = ifft(fft(subx, len) .* fft_h); % calc circle conv
+        y(i:i + L - 1) = y_k(M:M + L - 1); % save y
+    end
+    y = y(1:N + M - 1);
+end
+

大三上/信号处理原理/hw/matlab/2019/ex3/main.m

@@ -0,0 +1,96 @@
+% 读取数据
+[audio1, fs] = audioread('./data/1.wav');
+[audio2, fs] = audioread('./data/2.wav');
+[audio3, fs] = audioread('./data/3.wav');
+fprintf("read freq: %dHz\n", fs);
+total_len = 150000; % 取前150000个点
+audio1 = audio1(1 : total_len);
+audio2 = audio2(1 : total_len);
+audio3 = audio3(1 : total_len);
+
+% 原始的时域信号
+t = 0 : total_len - 1;
+figure(1);
+subplot(311);
+stem(t, audio1, '.');
+subplot(312);
+stem(t, audio2, '.');
+subplot(313);
+stem(t, audio3, '.');
+
+% 原始的频域信号
+figure(2);
+subplot(311);
+stem(t, abs(fft(audio1)), '.');
+subplot(312);
+stem(t, abs(fft(audio2)), '.');
+subplot(313);
+stem(t, abs(fft(audio3)), '.');
+
+% 1/3抽取和三倍上采样的下标
+sample_idx = 1 : 3 : total_len - 1;
+
+% 在时域进行1/3抽取
+audio1 = audio1(sample_idx);
+audio2 = audio2(sample_idx);
+audio3 = audio3(sample_idx);
+
+% 时域上进行上采样
+a1 = zeros(1, total_len);
+a2 = zeros(1, total_len);
+a3 = zeros(1, total_len);
+a1(sample_idx) = audio1;
+a2(sample_idx) = audio2;
+a3(sample_idx) = audio3;
+
+% 变换到频域
+f1 = fft(a1);
+f2 = fft(a2);
+f3 = fft(a3);
+
+% f是频域编码信号
+len = total_len / 6;
+f = zeros(1, total_len);
+f(1:len*3) 		 = [f1(1:len), 		 f2(1:len),		  f3(1:len)];
+f(len*3+1:len*6) = [f3(len+1:len*2), f2(len+1:len*2), f1(len+1:len*2)];
+
+% y是编码完成后的时域信号
+y = ifft(f);
+t = 0:total_len-1;
+figure(3);
+subplot(211);
+stem(t, abs(f), '.');
+subplot(212);
+stem(t, real(y), '.');
+
+% 频域解码
+f = fft(y);
+f1 = [f(1:len), 		f(len*5+1:len*6)];
+f2 = [f(len+1:len*2),	f(len*4+1:len*5)];
+f3 = [f(len*2+1:len*3), f(len*3+1:len*4)];
+t = 0:2*len-1;
+figure(4);
+subplot(311);
+stem(t, abs(f1), '.');
+subplot(312);
+stem(t, abs(f2), '.');
+subplot(313);
+stem(t, abs(f3), '.');
+
+% 恢复时域信号
+audio1 = real(ifft(f1));
+audio2 = real(ifft(f2));
+audio3 = real(ifft(f3));
+a1 = zeros(1, total_len);
+a2 = zeros(1, total_len);
+a3 = zeros(1, total_len);
+a1(sample_idx) = audio1;
+a2(sample_idx) = audio2;
+a3(sample_idx) = audio3;
+figure(5);
+subplot(311);
+stem(1:total_len, a1, '.');
+subplot(312);
+stem(1:total_len, a2, '.');
+subplot(313);
+stem(1:total_len, a3, '.');