Skip to content
Snippets Groups Projects
Select Git revision
  • c7bc2a98617d2578eca65f7921f9a4cb82112097
  • master default
  • dev
3 results

AssemblerMenuView.js

Blame
  • compressed_sensing.py 2.50 KiB
    import numpy as np
    import scipy
    import matplotlib.pyplot as plt
    
    
    def sample_two_sins(f1, f2, sample_times):
        sample_rads = 2 * np.pi * sample_times
        return np.sin(f1 * sample_rads) + np.sin(f2 * sample_rads)
    
    
    def compute_dct_matrix(n_samples):
        dct_matrix = np.zeros((n_samples, n_samples))
        root_one_over_n = np.sqrt(1.0 / n_samples)
        root_two_over_n = np.sqrt(2.0 / n_samples)
        for j in range(0, n_samples):
            dct_matrix[0, j] = root_one_over_n
        for i in range(1, n_samples):
            for j in range(0, n_samples):
                dct_matrix[i, j] = root_two_over_n * np.cos(np.pi * (2 * j + 1) * i / (2 * n_samples))
        return dct_matrix
    
    
    def compute_differences(recovered_dct, inverse_dct_matrix, sample_values):
        return sample_values - np.matmul(inverse_dct_matrix, recovered_dct)
    
    
    def loss_e(differences, subset_indices):
        return np.linalg.norm(differences[subset_indices])
    
    
    def grad_e(differences, dct_matrix, subset_indices):
        return -2 * np.matmul(dct_matrix[:,subset_indices], differences[subset_indices])
    
    
    if __name__ == "__main__":
        f1 = 697 # Hz
        f2 = 1209 # Hz
    
        # Part (a)
        sample_period = 0.01
        n_samples = 250
        sample_times = (sample_period / n_samples) * np.arange(n_samples)
        sample_values = sample_two_sins(f1, f2, sample_times)
        plt.plot(sample_times, sample_values)
        plt.savefig("fig_a.png")
        plt.close()
    
        # Part (b)
        dct_matrix = compute_dct_matrix(n_samples)
        dct = np.matmul(dct_matrix, sample_values)
        plt.plot(np.arange(n_samples), dct)
        plt.savefig("fig_b.png")
        plt.close()
    
        # Part (c)
        inverse_dct_matrix = np.transpose(dct_matrix)
        recovered_sample_values = np.matmul(inverse_dct_matrix, dct)
        plt.plot(sample_times, recovered_sample_values)
        plt.savefig("fig_c.png")
        plt.close()
    
        # Part (d)
        n_subsamples = 100
        np.random.seed(7250147)
        subset_indices = np.arange(n_samples)
        np.random.shuffle(subset_indices)
        subset_indices = subset_indices[:n_subsamples]
        subset_indices = np.sort(subset_indices)
        subset_sample_times = sample_times[subset_indices]
        subset_sample_values = sample_values[subset_indices]
        plt.plot(subset_sample_times, subset_sample_values)
        plt.savefig("fig_d.png")
        plt.close()
    
        # Part (e)
        recovered_dct = np.random.normal(0.0, 1.0, (n_samples))
        differences = compute_differences(recovered_dct, inverse_dct_matrix, sample_values)
        loss = loss_e(differences, subset_indices)
        grad = grad_e(differences, dct_matrix, subset_indices)
        print(grad)