Last active
March 15, 2017 09:28
-
-
Save FelipeCortez/432004cb91d0e955167d511da68c6997 to your computer and use it in GitHub Desktop.
FFT and reconstruction
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment