fourier.sci

function [H] = my_dft(signal)

    len = length(signal)
    H = zeros(1:len)

    for i = 1:len - 1
        for j = 1:len
            H(i) = H(i) + signal(j) * exp(-2 * %i * %pi * j * i / len)
        end
    end
endfunction

function [H] = my_fft(signal)

    disp('computing cooley-turkey')

    // cooley-turkey
    len = length(signal)

    disp(signal)
    disp(len)
    disp('^^^^^')

    if len == 1
        H = signal
        return
    end

    half_index = ceil(len / 2)


    even = zeros(1:half_index)
    odd = zeros(1:half_index)

    disp('computing odds and evens')

    even_k = 1
    odd_k = 1

    for i = 1:len

        if modulo(i, 2) == 0
           even(even_k) = signal(i)
           even_k = even_k + 1
        else
           odd(odd_k) = signal(i) 
           odd_k = odd_k + 1
        end
    end

    output = zeros(1:len)

    disp('Starting to define submatrix')

    output(1:half_index)       = my_fft(even)
    output(half_index + 1:len) = my_fft(odd)

    disp('Computing the final result')

    for k = 1:half_index
        t = output(k)                
        exponent = exp(-2 * %i * %pi * k / len)

        output(k) = t + exponent * output(k + half_index)
        output(k + half_index) = t - exponent * output(k + half_index)
    end

    H = output
    disp('terminating')    
endfunction

function [H] = compute_mse(x1, x2)
    H = (x1 - x2) ** 2 / length(x1)
endfunction

// result = my_fft(ones(1:16))

signal = sin(2 * %pi * linspace(0, 10, 128))

mine = abs(my_dft(signal))
theirs = abs(fft(signal))


mine = cat(2, [0], mine)(1:length(mine))
freq = (0:length(mine) - 1) * length(signal) / 10 / length(mine)

ball_size = 10
red = [144, 0, 0]
green = [0, 176, 0]


subplot(211) ; freq_len = round(length(freq) / 2)
plot2d(freq, mine, rect=[0, -10, round(max(freq) / 2), max(mine) * 1.5], style=color('red')) ; title("My fft")
scatter(freq(1:freq_len), mine(1:freq_len), ball_size)

// subplot(412)
// plot2d(freq, theirs, rect=[0, -10, round(max(freq) / 2), max(theirs) * 1.5], style=color('green')) ; title("Canonical dft")
// scatter(freq(1:freq_len), theirs(1:freq_len), ball_size)

subplot(212)
mse = (theirs - mine) ** 2 / length(theirs)
plot2d(freq, mse, rect=[0, min(mse) - 2.5, round(max(freq) / 2), max(mse) + 2.5], style=color('green')) ; title("Mean Square Error")
scatter(freq(1:freq_len), mse(1:freq_len), ball_size)


// subplot(414)
// plot2d(signal, rect=[0, min(signal) * 1.5, length(signal), max(signal) * 1.5]) ; title("Original signal")
// scatter(1:length(signal), signal, ball_size * 0.5)
// 

task2.sci

A = 1

function [H] = sin1(x)
    H = A * sin(2 / 3 * 2 * %pi * x)  // T = 1.5      // red
endfunction

function [H] = sin2(x)
    H = A * sin(2.5 * 2 * %pi * x)    // T = 4        // green
endfunction

function [H] = sin3(x)
    H = A * sin(5 * 2 * %pi * x)      // T = 0.2      // blue // T = 10 / 3 f = 0.3 would work tho
endfunction



function [H] = compute(x)  
    H = sin1(x) + sin2(x) + sin3(x) // + sin(1 * 2 * %pi * x) // sin(0.1 * 2 * %pi * x) 
endfunction


if ~exists("import_definitions", "local") then
    import_definitions = 1 == 0
end

if import_definitions
    disp('Definitions are imported')
    return
end

interval = 10.5 // time period
orig_Fs = 40 // .045  // > 3.5 * 2
orig_samples = orig_Fs * interval

ball_size = 10


x = linspace(0, interval, orig_samples); // (0:255); // 
func = compute(x);


std = 411
sanity = std // std // 111

Xs = linspace(0, interval, interval * 100);
Ys = compute(Xs);



subplot(sanity); plot2d(Xs, Ys, rect=[0, -4, 10.5, 4]); title('orig signal');
// subplot(sanity); scatter(Xs, Ys, ball_size * 0.1);


subplot(412); plot2d(Xs, sin1(Xs), rect=[0, -1.5, 10.5, 1.5], style=color('red')); title('constituent frequencies');
subplot(412); plot2d(Xs, sin2(Xs), style=color('green')); 
subplot(412); plot2d(Xs, sin3(Xs), style=color('blue')); 

if sanity == 111
    return
end

subplot(413); scatter(x, func, ball_size);
subplot(413); plot2d(x, func, rect=[0, -4, 10.5, 4]); title('sampled signal');


fourier_image_orig = floor(abs(fft(func)));
orig_freq = (0:length(fourier_image_orig) - 1) * orig_samples / interval / length(fourier_image_orig);

fourier_image_orig = fourier_image_orig(1:round(length(fourier_image_orig) / 2))
orig_freq = orig_freq(1:round(length(orig_freq) / 2))

miny_orig = max(fourier_image_orig) + 0.3 * max(fourier_image_orig)
minx_orig = max(orig_freq) / 2 


scf(0);
subplot(414)

scatter(orig_freq, fourier_image_orig, ball_size);
plot2d(orig_freq, fourier_image_orig, rect=[0, 0, minx_orig, miny_orig]); title('canonical fft');

// subplot(222);
// plot2d(orig_freq, abs(my_fft(func))); title('my fft');

// to show these values in the terminal

interval = interval
orig_samples = orig_samples

task3.sci

import_definitions = 1 == 1 ; exec('task2.sci')


interval = 20
Fs = 20
samples_n = Fs * interval

x = linspace(0, interval, samples_n); // (0:255); // 

func = compute(x)
func_padded = resize_matrix(func, 1, samples_n * 2)
fourier_padded = abs(fft(func_padded))

fourier = abs(fft(func))

freq = (0:length(fourier) - 1) * samples_n / interval / length(fourier);
freq_padded = (0:length(fourier_padded) - 1) * samples_n / interval / length(fourier_padded);




miny_orig = max(fourier) + 0.3 * max(fourier)
minx_orig = max(freq) / 2

subplot(211)
plot2d(freq_padded, fourier_padded, rect=[0, 0, minx_orig, miny_orig]) ; title("Frequency of the padded sampled signal")
subplot(212)
plot2d(freq, fourier, rect=[0, 0, minx_orig, miny_orig]) ; title("Frequency of the orig sampled signal")

task4.sci

function [H] = cosine_signal(x, A, f)
    H = A * cos(f * 2 * %pi * x);
endfunction


interval = 1 // seconds
Fs = 200     // S/s // 570

samples_n = Fs * interval
x = linspace(0, interval, samples_n)

signal1 = cosine_signal(x, 0.5, 190);
signal2 = cosine_signal(x, 2, 10);


freq1 = (0:length(signal1) - 1) * Fs / length(signal1);
freq2 = freq1;

fourier1 = abs(fft(signal1));
fourier2 = abs(fft(signal2));


miny1 = max(fourier1) + 0.3 * max(fourier1)
minx1 = max(freq1) / 2

miny2 = max(fourier2) + 0.3 * max(fourier2)
minx2 = max(freq2) / 2

subplot(411)
plot2d(x, signal1, style=color('magenta'))

orig_x = linspace(0, interval, 8196)
plot2d(orig_x, cosine_signal(orig_x, 0.5, 190) , style=color('red'))
subplot(412)
plot2d(x, signal2, style=color('green'))

subplot(413)
plot2d(freq1, fourier1, rect=[0, 0, minx1, miny1]) ; title("Frequency of sinusoid 1")
subplot(414)
plot2d(freq2, fourier2, rect=[0, 0, minx2, miny2]) ; title("Frequency of sinusoid 2")

samples_n = samples_n