clear;                  % Clear Workspace and Command Window
clc;

UseTwosComplement = 1;
StraighBinary = 0;
fs = 40 * 10^6 / 4;
NumberOfADCBits = 14;
Fullscale = 1;              % +- 1V
FixEndian = 1;

    %open file and read data
    fid = fopen('Output.hex');
       FileData = fread(fid,'uint16=>uint16'); 
    fclose(fid);


if FixEndian == 1
    FileData = swapbytes(FileData);
end;
 
FileData = bitand(FileData, uint16(hex2dec('3FFF')));       % take care of overflow
    
U_t = cast (FileData, 'double');

if UseTwosComplement == 1
    U_t = TwosComplementToDecimal(U_t, NumberOfADCBits);
end;

if StraighBinary == 1
    U_t = U_t + 2^13-1;
else
    %U_t = U_t / 2^(NumberOfADCBits-1)*Fullscale/30;
end;

t_A = [0:(1/fs):(length(U_t)-1)*(1/fs)];
figure(3); plot(t_A, U_t, 'blue'); grid on;

% FFT / DFT
fn = fs / 2;                                % Nyquistfrequenz
Nfft = length(U_t);
df = fs / Nfft;                             % Correct frequency spacing

f_axis = (-Nfft/2 : Nfft/2-1) * df;

% FFT
SequenceSpectrum = fft(U_t, Nfft)/(Nfft/2);     % Get FFT and Normalize
SequenceSpectrum = fftshift(SequenceSpectrum);

% FFT Angle
SequenceSpectrum_Angle = unwrap(angle(SequenceSpectrum));
SequenceSpectrum_Angle = SequenceSpectrum_Angle * 180 / pi;

% Remove Noise
%SequenceSpectrum((SequenceSpectrum<100)) = 0;

figure(1);
subplot(3,1,1);  stem(f_axis, real(SequenceSpectrum)); title('Real part'); xlabel('f / Hz'); ylabel(''); grid on;
subplot(3,1,2);  stem(f_axis, imag(SequenceSpectrum)); title('Imagenary part');  xlabel('f / Hz'); ylabel('');grid on;
subplot(3,1,3);  stem(f_axis, abs(SequenceSpectrum)); title('Magnitude'); xlabel('f / Hz'); ylabel('');  grid on;


Y = ifftshift(SequenceSpectrum);
y = ifft((fftshift(Y)), Nfft);

t = [0:(1/fs):(length(real(y))-1)*(1/fs)];

figure(2);
plot(t, real(y)); title('Synthesis'); xlabel('t / s'); ylabel('');  grid on;

fprintf('NFFT %i\n', Nfft);
fprintf('df %i Hz\n', df);