FFT and reconstruction

lab3.sce

function [magnitude, phase] = fourierSeries(x)
    // FOURIERSERIES - Return the magnitude and phase of each
    // sinusoidal component in the Fourier series expansion for
    // the argument, which is interpreted as one cycle of a
    // periodic signal.  The argument is assumed to have an
    // even number p of samples. The first returned value is an
    // array containing the amplitudes of the sinusoidal
    // components in the Fourier series expansion, with
    // frequencies 0, 1/p, 2/p, ... 1/2. The second returned
    // value is an array of phases for the sinusoidal
    // components. Both returned values are arrays with length
    // (p/2)+1.
    p = length(x);
    f = fft(x)/p;
    magnitude(1) = abs(f(1));
    upper = p/2;
    magnitude(2:upper) = 2*abs(f(2:upper));
    magnitude(upper+1) = abs(f(upper+1));
    phase(1) = atan(f(1));
    phase(2:upper) = atan(f(2:upper));
    phase(upper+1) = atan(f(upper+1));
endfunction

function x = reconstruct(magnitude, phase)
    // RECONSTRUCT - Given a vector of magnitudes and a vector
    // of phases, construct a signal that has these magnitudes
    // and phases as its discrete Fourier series coefficients.
    // The arguments are assumed to have odd length, p/2 + 1,
    // and the returned vector will have length p.
    t = [0 : 1/rate : 1 - 1/rate]
    p = (length(magnitude) - 1) * 2
    x = zeros(p) + magnitude(1)

    for i = 2:length(magnitude)
        x = x + (magnitude(i) * cos(i * 2 * %pi * t + atan(imag(phase(i)), real(phase(i)))))
    end
endfunction

rate = 8000
t = [0 : 1/rate : 1 - 1/rate]
x = sin(2 * %pi * 800 * t)

// 1) --------

function [] = ex1()    clf()
    [A, phi] = fourierSeries(x)
    p = length(x)
    freqs = [0 : rate/p : rate/2]
    plot2d3(freqs, A)
endfunction

// 2) --------

function [] = ex2()
    clf()
    y = sin(2 * %pi * 800 * (t .* t));
    //y = sin(2 * %pi * 800 * t);
    [A, phi] = fourierSeries(y)
    subplot(211)
    plot2d3(t, y)
    //plot2d3(freqs, A)
    subplot(212)
    rec = reconstruct(A, phi)
    plot2d3(rec)
endfunction

// 3) --------

function [] = ex3()
    clf()
    wave_1 = 1/2 * cos(2 * %pi * 790 * t)
    wave_2 = 1/2 * cos(2 * %pi * 810 * t)
    wave_sum = wave_1 + wave_2
    plot2d3(t(1:800), wave_sum(1:800))
endfunction

// 4)

function [] = ex4()
    clf()
    [A, phi] = fourierSeries(wave_sum)
    p = length(wave_sum)
    freqs = [0 : rate/p : rate/2]
    plot2d3(freqs, A)
endfunction