diff --git a/matlab/DFT_time2freq.m b/matlab/DFT_time2freq.m
new file mode 100644
index 0000000..79baac1
--- /dev/null
+++ b/matlab/DFT_time2freq.m
@@ -0,0 +1,16 @@
+function f_val = DFT_time2freq( t, val, freq )
+% val = FFT_time2freq( t, val, freq )
+%
+% computes the DFT at the given frequencies
+
+if numel(t) ~= numel(val)
+    error 'numel(t) ~= numel(val)'
+end
+
+f_val = zeros(1,numel(freq));
+for f_idx=1:numel(freq)
+    f_val(f_idx) = sum( val .* exp( -1i * 2*pi*freq(f_idx) * t ) );
+end
+f_val = f_val / numel(t);
+
+f_val = f_val * 2; % single-sided spectrum
diff --git a/matlab/FFT_time2freq.m b/matlab/FFT_time2freq.m
index cdfad45..ece7fd9 100644
--- a/matlab/FFT_time2freq.m
+++ b/matlab/FFT_time2freq.m
@@ -1,9 +1,10 @@
 function [f,val] = FFT_time2freq( t, val )
 
-dt=t(2)-t(1);
-val = [val zeros(1,5000)];
-L=numel(val);
-f = (0:L-1)/L/dt;
-f = f(1:floor(L/2));
-val = 2*fft(val)/L;
-val = val(1:floor(L/2));
+dt=t(2)-t(1); % timestep
+L=numel(val); % signal length
+NFFT = 2^nextpow2(L); % next power of 2 (makes fft fast)
+%very fine freq resolution... NFFT = NFFT+100000;
+val = fft( val, NFFT)/L;
+f = 1/(2*dt) * linspace(0,1,NFFT/2+1);
+
+val = 2*val; % single-sided spectrum