Determine original function after fourier transform

Dear readers,

I am currently working on a school assignment where i get a function, 20000 samples long and i have to determine the original frequences in this signal. There is also noise in the signal. To do this i first fourier transformed the signal, and created the phase and amplitude spectra. By setting a minimum in the amplitude spectra i hope to reduce the noise. Now i want to return to a signal in the time domain, since i know what the amplitude and phase of each frequence is, this should be possible. To do this i tried to use the symsum command, but it complains about the matrix dimentions.

I have added my current code below, note that all the remarks are in dutch. The last part of the code is where i run into problems. P1 and fase are both arrays with the same length, ~16000 samples currently.

note: i am not a native speaker, feel free to ask any questions.

clear all

close all

clc

%Definieren van samplefrequentie(Fs),

%aantal nullen(2^n), lengte van het

%uiteindelijke signaal(Lsig), het

%window(W) dat toegepast wordt en

%de uiteindelijke array(signaal)

%n moet minimaal 15 zijn

Fs = 5000;

n = 15;

Lsig = 2^n;

W = window(@chebwin,20000);

signaal = zeros(Lsig,1);

%inladen van het signaal

load signaal

%sig = sig.*W;

%Het signaal plaatsen in de array

%gevuld met nullen

for k = 1:20000

signaal(k) = sig(k);

end

%vermenigvuldigen van het signaal

%met de window

%signaal = signaal.*W;

%Fourier transformatie toepassen

%op het signaal

Y = fft(signaal,Lsig);

%fase van

fase1 = angle(Y(1:Lsig/2+1));

fase = fase1./pi();

%amplituderespons van het signaal

%bepalen en halveren omdat het

%signaal zich daarna herhaalt

P2 = abs(Y/Lsig);

P1 = 20.*P2(1:Lsig/2+1);

%de amplituderespons die lager zijn dan de

%gegeven waarde bij het if statement gelijk

%stellen aan 0

for z = 1:length(P1)

if P1(z) < 1.8*10^-2

P1(z) = 0;

else

end

end

%x-as definieren en plots maken

f = Fs*(0:(Lsig/2))/Lsig;

figure(1)

semilogy(f,P1)

xlabel('f (Hz)')

ylabel('|P1(f)| (dB)')

title('Amplituderespons als functie van de frequentie')

figure(2)

plot(f,fase)

xlabel('f (Hz)')

ylabel('faseverschil')

%oorspronkelijke functie terug halen

x = 1:20000;

sig0 = zeros(20000,1);

for l = 1:length(P1);

sig0(x) = symsum(P1(l).*sin(2.*pi().*f.*x + fase(l).*x), x, 1, 20000);

end