🚌 programs, results

ImTest.m

@@ -0,0 +1,3 @@
+A = imread('crysis.jpg');
+A = double(A) / 255;
+imshow(A);

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2020 nitrece
+Copyright (c) 2012-20 Subhajit Sahu
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

README.md

@@ -1,2 +1,7 @@
-# digital-signal-processing-laboratory
-Digital signal processing (DSP) is the use of computers or specialized hardware, to perform a wide variety of signal processing operations.
+Digital signal processing (DSP) is the use of computers or specialized
+hardware, to perform a wide variety of signal processing operations.
+
+**Course**: Digital Signal Processing Laboratory, [Monsoon 2012]<br>
+**Taught by**: Prof. Sumit Saha
+
+[Monsoon 2012]: https://github.com/nitrece/semester-5

circonv.m

@@ -0,0 +1,15 @@
+function y = circonv(x, h)
+xSz = length(x);
+hSz = length(h);
+ySz = max([xSz hSz]);
+x = [x zeros(1, ySz - xSz)];
+h = [h zeros(1, ySz - hSz)];
+
+% Generate Xc matrix
+Xc = zeros(ySz, ySz);
+for i = 0 : (ySz - 1)
+    Xc(1 + i, :) = [x(1, (ySz + 1 - i) : ySz) x(1, 1 : (ySz - i))];
+end
+
+y = (Xc' * h')';
+end

conv_overlap_add.m

@@ -0,0 +1,28 @@
+function y = conv_overlap_add(x, h)
+L = 4;
+M = 4;
+N = L + M - 1;
+xLen = length(x);
+hLen = length(h);
+T = ceil(xLen / L) + 1;
+H = [h zeros(1, N - hLen)];
+
+% Perform Overlap-save method
+X = zeros(T, N);
+Y = X;
+for i = 0 : (T-1)
+    if(i == T-1)
+        X(1+i, :) = [x(1, (1+i*L):xLen) zeros(1, T*L - xLen) zeros(1, M-1)];
+    else
+        X(1+i, :) = [x(1, (1+i*L):(i*L+L)) zeros(1, M-1)];
+    end
+Y(1+i, :) = circonv(X(1+i, :), H);
+end
+
+% Merge useless output and generate output
+y = zeros(1, T*L+M);
+for i = 0 : (T-1)
+    y(1, (1+i*L):(i*L+N)) =  y(1, (1+i*L):(i*L+N)) + Y(1+i, 1:N);
+end
+y = y(1, 1 : (xLen + hLen - 1));
+end

conv_overlap_save.m

@@ -0,0 +1,30 @@
+function y = conv_overlap_save(x, h)
+L = 4;
+M = 4;
+N = L + M - 1;
+xLen = length(x);
+hLen = length(h);
+T = ceil(xLen / L) + 1;
+H = [h zeros(1, N - hLen)];
+
+% Perform Overlap-save method
+X = zeros(T, N);
+Y = X;
+for i = 0 : (T-1)
+    if(i == 0)
+        X(1, :) = [zeros(1, M-1) x(1, 1:L)];
+    elseif(i == T-1)
+        X(1+i, :) = [X(i, (1+L):N) x(1, (1+i*L):xLen) zeros(1, T*L - xLen)];
+    else
+        X(1+i, :) = [X(i, (1+L):N) x(1, (1+i*L):(i*L+L))];
+    end
+Y(1 + i, :) = circonv(X(1+i, :), H);
+end
+
+% Discard useless output and generate output
+y = zeros(1, T*L);
+for i = 0 : (T-1)
+    y(1, (1+i*L):(i*L+L)) = Y(1+i, M:end);
+end
+y = y(1, 1 : (xLen + hLen - 1));
+end

convolve.m

@@ -0,0 +1,9 @@
+function y = convolve(x, h)
+xSz = length(x);
+hSz = length(h);
+ySz = xSz + hSz - 1;
+y = zeros(1, ySz);
+for i = 1 : xSz
+    y(1, i:(i+hSz-1)) = y(1, i:(i+hSz-1)) + x(1, i) * h;
+end
+end

dft_circonv.m

@@ -0,0 +1,9 @@
+function y = dft_circonv(x, h)
+% Circular Convolution using DFT
+xLen = length(x);
+hLen = length(h);
+yLen = max([xLen hLen]);
+x = [x zeros(1, yLen - xLen)];
+h = [h zeros(1, yLen - hLen)];
+y = sig_idft(sig_dft(x) .* sig_dft(h));
+end

dft_conv_test.m

@@ -0,0 +1,24 @@
+clc;
+clear all;
+close all;
+
+t0 = linspace(0, 6*pi, 1000);
+t1 = linspace(0, 2*pi, 100);
+x = cos(t0);
+h = cos(t1);
+y0 = real(dft_convolve(x, h));
+subplot(2,2,1);
+plot(x, 'LineWidth', 2);
+title('Time Domain Signal x(n)');
+xlabel('Sample Number');
+ylabel('Amplitude');
+subplot(2,2,2);
+plot(h, 'LineWidth', 2);
+title('Impulse Response h(n)');
+xlabel('Sample Number');
+ylabel('Amplitude');
+subplot(2,2,3:4);
+plot(y0, 'LineWidth', 2);
+title('Ouput Signal y(n) using DFT and IDFT');
+xlabel('Sample Number');
+ylabel('Amplitude');

dft_convolve.m

@@ -0,0 +1,9 @@
+function y = dft_convolve(x, h)
+% Linear Convolution using DFT
+xLen = length(x);
+hLen = length(h);
+yLen = xLen + hLen - 1;
+x = [x zeros(1, yLen - xLen)];
+h = [h zeros(1, yLen - hLen)];
+y = sig_idft(sig_dft(x) .* sig_dft(h));
+end

dft_test.m

@@ -0,0 +1,18 @@
+clc;
+clear all;
+close all;
+
+t0 = linspace(0, 1, 1000);
+x = cos(2*pi*30*t0);
+y = sig_dft(x);
+
+subplot(2,1,1);
+plot(t0, x, 'LineWidth', 2);
+title('Time Domain Signal');
+xlabel('Time');
+ylabel('Amplitude');
+subplot(2,1,2);
+plot(linspace(-500, 500, 1000), fftshift(abs(y)), 'LineWidth', 2);
+title('Frequency Domain Signal');
+xlabel('Frequency');
+ylabel('Amplitude');

discretize.m

@@ -0,0 +1,5 @@
+function y = discretize(x, amp, bits)
+m_bits = bits - 1;
+y = round((x / amp) * (2 ^ m_bits - 1));
+y = (y / (2 ^ m_bits - 1)) * amp;
+end

fft_convolve.m

@@ -0,0 +1,8 @@
+function y = fft_convolve(x, h)
+xLen = length(x);
+hLen = length(h);
+yLen = xLen + hLen - 1;
+x = [x zeros(1, yLen - xLen)];
+h = [h zeros(1, yLen - hLen)];
+y = ifft(fft(x) .* fft(h));
+end

fir_convolve.m

@@ -0,0 +1,12 @@
+function y = fir_convolve(x, h)
+lenX = length(x);
+lenH = length(h);
+lenY = lenX + lenH - 1;
+x = [x zeros(1, lenH-1)];
+y = zeros(1, lenY);
+z = zeros(1, lenH);
+for i = 1 : lenY
+    z = [x(1, i) z(1, 1:(lenH-1))];
+    y(1, i) = sum(z .* h);
+end
+end

quantization_err_in_conv.m

@@ -0,0 +1,22 @@
+clc;
+clear all;
+close all;
+
+t0 = linspace(0, 6*pi, 1000);
+t1 = linspace(0, 2*pi, 100);
+x = cos(t0);
+h = cos(t1);
+y0 = conv(x, h);
+
+pLen = 128;
+p = zeros(1, pLen);
+for i = 2 : pLen
+    y1 = quantized_convolve(x, 1, h, 1, i);
+    p(1, i) = sum((y1 - y0) .^ 2) / sum(y0 .^ 2);
+end
+p(1, 1) = p(1, 2);
+
+plot(10 * log10(p));
+title('Quantization Error in Convolution');
+xlabel('Bits used');
+ylabel('dB Error');

quantized_convolve.m

@@ -0,0 +1,16 @@
+function y = quantized_convolve(x, x_amp, h, h_amp, bits)
+lenX = length(x);
+lenH = length(h);
+lenY = lenX + lenH - 1;
+z_amp = x_amp * h_amp;
+y_amp = z_amp * lenH;
+x = discretize(x, x_amp, bits);
+h = discretize(h, h_amp, bits);
+x = [x zeros(1, lenH-1)];
+y = zeros(1, lenY);
+z = zeros(1, lenH);
+for i = 1 : lenY
+    z = [x(1, i) z(1, 1:(lenH-1))];
+    y(1, i) = discretize(sum(discretize(z .* h, z_amp, bits)), y_amp, bits);
+end
+end

sig_dft.m

@@ -0,0 +1,10 @@
+function y = sig_dft(x)
+% Discrete Fourier Transform
+N = length(x);
+y = zeros(1, N);
+for k = 0 : (N-1)
+    for n = 0 : (N-1)
+        y(1, 1 + k) = y(1, 1 + k) + x(1, 1 + n) * exp(-1j * (2*pi/N) * n * k);
+    end
+end
+end

sig_idft.m

@@ -0,0 +1,11 @@
+function y = sig_idft(x)
+% Inverse Discrete Fourier Transform
+N = length(x);
+y = zeros(1, N);
+for n = 0 : (N-1)
+    for k = 0 : (N-1)
+        y(1, 1 + n) = y(1, 1 + n) + x(1, 1 + k) * exp(1j * (2*pi/N) * n * k);
+    end
+    y(1, 1 + n) = y(1, 1 + n) / N;
+end
+end

test_circonv.m

@@ -0,0 +1,25 @@
+% Circular Convolution
+
+% Calculate Values and perform Circular Convolution
+t = linspace(0, 1, 1024);
+x = sin(2*pi*16*t);
+h = log(abs(cos(2*pi*8*(t-0.5))./(pi*8*(t-0.5)))) / 5;
+y = circonv(x, h);
+
+% Plot Figures
+figure;
+subplot(2, 2, 1);
+plot(t, x, 'LineWidth', 2);
+title('Input Signal x(n)');
+xlabel('Time');
+ylabel('Amplitude');
+subplot(2, 2, 2);
+plot(t, h, 'LineWidth', 2);
+title('Impulse Response h(n)');
+xlabel('Time');
+ylabel('Amplitude');
+subplot(2, 2, 3:4);
+plot(t, y, 'LineWidth', 2);
+title('Circular-Convolved Signal y(n) = x(n) * h(n)');
+xlabel('Time');
+ylabel('Amplitude');

test_dft_circonv.m

@@ -0,0 +1,25 @@
+% Circular Convolution using DFT
+
+% Calculate Values and perform Circular Convolution
+t = linspace(0, 1, 1024);
+x = sin(2*pi*16*t);
+h = sin(2*pi*8*(t-0.5))./(pi*8*(t-0.5));
+y = dft_circonv(x, h);
+
+% Plot Figures
+figure;
+subplot(2, 2, 1);
+plot(t, x, 'LineWidth', 2);
+title('Input Signal x(n)');
+xlabel('Time');
+ylabel('Amplitude');
+subplot(2, 2, 2);
+plot(t, h, 'LineWidth', 2);
+title('Impulse Response h(n)');
+xlabel('Time');
+ylabel('Amplitude');
+subplot(2, 2, 3:4);
+plot(t, y, 'LineWidth', 2);
+title('Circular-Convolved Signal y(n) = x(n) * h(n)');
+xlabel('Time');
+ylabel('Amplitude');

using_convolution.m

@@ -0,0 +1,22 @@
+%
+% Using Convolution
+%
+clc;
+clear all;
+close all;
+
+t = linspace(0, 1, 1000);
+x = sin(6*pi*t);
+h = mod(t, 0.1);
+y = conv(x, h);
+
+figure;
+subplot(5, 2, [1 3]);
+plot(x);
+title('Input Signal (x)');
+subplot(5, 2, [7 9]);
+plot(h);
+title('Impulse Response (h)');
+subplot(5, 2, [4 6 8]);
+plot(y);
+title('Output Signal (y)');