[data] = loadQucsDataSet("fft_test.dat")
showQucsDataSet(data)
[Vout] = getQucsVariable(data,"out.Vt");
[time]  = getQucsVariable(data,"time");
Ts = time(2) - time(1);
fs = 1/Ts;
N = length(Vout);%Ensure that N >> fs/2.
y_fft = fft(Vout, N);%fft returns a vector pi rad shifted. Since we're dealing with real signals, the spectrum is symmetric around f=0. Let's cut the first half...
y_fft = y_fft(1 : round(N/2));%Just cut the first half. No fftshift is needed.
y_fft = 2*y_fft/N;%Now the result must be normalized wrt the number of points. OTOH, since we've just removed the second half of the spectrum, we must multiply the result by 2.
freq = linspace(0, fs/2, N/2);%Prepare the freq axis...
clf();
y_abs = abs(y_fft);

%Check the result is correct
[max_val, posx] = max(y_abs);
max_val
f_tone = freq(posx);
f_tone

%Plot the result
span = 3e3;%3kHz
fc = 1e3;
plot(freq, y_abs)
title("Single-sided spectrum");
xlabel("Frequency (Hz)");
ylabel("Amplitude (V)");
grid on
axis ([fc-span/2 fc+span/2])