What is wrong with my FFT code?

I have this code, the plot looks correct but the values don't come out right. I even tested an ideal sine wave, but it's the same. which means the code is wrong, but I can't figure out which part of the code is wrong. it looks correct to me. please help me out;

% Parameters

N = 8192; % Number of samples for FFT

fs = 10e9; % Sampling frequency (10 GHz)

harmonics = 3; % Number of harmonics to consider

% Extract 8192 samples from OUTPUT (assuming OUTPUT is already defined)

x = OUTPUT(1:N);

% Normalize data to unit magnitude

x = x / max(abs(x));

% Perform FFT

X = fft(x);

% Frequency axis

f_axis = (0:N-1) * fs / N; % Frequency axis from 0 to fs

% Remove DC component

X(1) = 0;

% Find dominant frequency component

[~, max_idx] = max(abs(X)); % Find the index of the maximum magnitude

% Calculate estimated signal frequency

f_signal = f_axis(max_idx); % Frequency corresponding to the maximum magnitude

% Calculate SNR correctly

signal_power = abs(X(max_idx))^2; % Power of the signal tone

noise_power = sum(abs(X).^2) - signal_power; % Total noise power (excluding signal tone)

SNR = 10*log10(signal_power / noise_power); % Signal-to-Noise Ratio in dB

% Calculate ENOB

ENOB = (SNR - 1.76) / 6.02;

% Calculate SFDR

SFDR_power = max(abs(X))^2; % Power of the largest spurious component (excluding DC)

SFDR = 10*log10(SFDR_power / signal_power); % Spurious Free Dynamic Range in dB

% Calculate THD correctly

THD_power = sum(abs(X(2:end)).^2); % Total harmonic distortion power (excluding signal and DC)

THD = 10*log10(THD_power / signal_power); % Total harmonic distortion in dB

% Plot FFT magnitude spectrum (folded, centered around fs/2)

figure;

plot(f_axis / 1e6, 20*log10(abs(X) / N)); % Plot in MHz and dB

xlabel('Frequency (MHz)');

ylabel('Magnitude (dB/Hz)');

title('FFT Magnitude Spectrum');

grid on;

xlim([0 fs/2/1e6]); % Limit x-axis to 0 to fs/2 in MHz

% Mark signal tone and harmonics on the plot

hold on;

plot(f_signal / 1e6, 20*log10(abs(X(max_idx)) / N), 'ro', 'MarkerSize', 10, 'LineWidth', 2); % Mark signal tone

for i = 1:harmonics

harmonic_freq = i * f_signal;

[~, harmonic_bin] = min(abs(f_axis - harmonic_freq)); % Find nearest bin index

plot(f_axis(harmonic_bin) / 1e6, 20*log10(abs(X(harmonic_bin)) / N), 'go', 'MarkerSize', 10, 'LineWidth', 2); % Mark harmonics

end

hold off;

% Create legend for ENOB, SNR, SFDR, THD without identifiers

legend_str = {['ENOB: ', num2str(ENOB), ' bits'], ['SNR: ', num2str(SNR), ' dB'], ...

['SFDR: ', num2str(SFDR), ' dB'], ['THD: ', num2str(THD), ' dB']};

legend(legend_str, 'Location', 'best', 'AutoUpdate', 'off'); % AutoUpdate off to prevent adding previous legend entries

% Display results

disp(['Estimated Signal Tone Frequency: ', num2str(f_signal / 1e6), ' MHz']);

disp(['ENOB: ', num2str(ENOB), ' bits']);

disp(['SNR: ', num2str(SNR), ' dB']);

disp(['SFDR: ', num2str(SFDR), ' dB']);

disp(['THD: ', num2str(THD), ' dB']);

[Paul]

Hi Ehsan,

What do you think is incorrect and why?

[Ehsan]

One of the most likely culprits is the signal power, if the signal power is not calculated properly, then everything else is wrong, as they involve signal power.

Another possiblity is that the values are not being calculated with correct units. But I just can't figure out why it could be the case.

Answer 1

[Umar]

Hi Ehsan,

You asked, I can't figure out which part of the code is wrong. it looks correct to me. please help me out;

Upon reviewing the code, I have identified a potential bug that could be causing the incorrect values. The bug lies in the calculation of the noise power, which is used to calculate the SNR. The line noise_power = sum(abs(X).^2) - signal_power; is incorrect because it subtracts the signal power from the total power, assuming that the remaining power is noise. However, this assumption is incorrect, as it includes the power of the harmonics and other frequency components. So, Instead of subtracting the signal power from the total power, we will subtract the power of the signal tone and the power of the harmonics from the total power. This will give us the correct noise power, which can then be used to calculate the SNR accurately.

Replace the line noise_power = sum(abs(X).^2) - signal_power; with the following code:

% Calculate noise power correctly noise_power = sum(abs(X).^2) - signal_power - sum(abs(X(harmonic_bin)).^2);

By subtracting the power of the harmonics (sum(abs(X(harmonic_bin)).^2)) from the total power, we exclude their contribution from the noise power calculation.

Hope this will help resolve your problem.

[Ehsan]

No, this answer does not seems to be correct. also logicaly speaking if we remove the signal tone from the total power spectral density, what is left is the harmonics and noise, so it makes sense to calculate noise that way. either way, I tried that change of code, and it did not work.

Answer 2

[Umar]

Hi Eshan,

The proposed correction to the noise power calculation by subtracting the power of the signal tone and harmonics from the total power seems logical. However, you mentioned that implementing this change did not resolve the issue. It's essential to ensure that the harmonic_bin variable correctly identifies the harmonics' frequency components. Double-check the indexing and calculation of the harmonics' power to accurately exclude their contribution from the noise power calculation. Additionally, verify that the signal_power calculation is correct and consistent with the signal tone's power. Debugging the code step by step and inspecting intermediate results can help pinpoint the exact source of the discrepancy.