The power of a sinusoidal component in noisy PSD

Question

I am very new to signal processing and I have been doing a lot of reading on it but I don't seem to understand everything. I'm hoping you can help my confusion a little.

I am simulating a lidar signal that has a lot of noise and the end goal is to get the power of the sinusoidal component. My algorithm for getting the spectrum that I am simulating is as follows:

1. capture signal window
2. apply a window function to it
3. do zero padding on a copy of the windowed signal
4. calculate the FFT of windowed signal and the zero padded signal
5. Calculate the norm squared of the two FFT outputs
6. Sum 10000 of those to degrease noise in spectrums
7. Multiply the spectrum by two to get a one sided FFT. divide by 10000 to get average. divide by window length squared (this I don't fully understand) and divide by window function constant (explained below)

My simulation matlab code is below.

My first problem is figuring out the power of a signal in low noise signal without window function. I understand the power of a sinusoidal signal is (amplitude^2)/2. So in this case the amplitude is 0.5 and power is 0.125. When I zero padd the signal the peak amplitude of my spectrum is very close to 0.125 but on the non zero padded signal the peak is not there. This could be due to the signal power being spread out on different bins but I jus can't figure out how to calculate the power from that. Here is the figure of freq1:

My second question concerns the noisier environment. In that situation the zero padded spectrum shows the amplitude as 0.125 as well but the non zero padded signal is raised way above the 0.125 and thus the amplitude of the peak is much higher. I could somehow normalize the spectrum but I don't know how exactly and what is the logical explanation of this process. Picture of this high noise signal at the same frequency peak is here:

My third question involves the window function. When I apply windowing to the signal it's amplitude in the fft magnitude spectrum is scaled by some constant (for example hanning windowed spectrum amplitude is 0.5 of not windowed spectrum). So to get the magnitude spectrum accurate I would divide by this constant. How does this work in the power spectrum? Do I just square the constant?

Fourth question. Why is the peak of my non zero padded signal no on the frequency bin it should be? Am I doing something wrong on the non zero padded signal analysis?

These all combined. How would I calculate the power of the sinusoidal component of freq1 in the case where there is a lot of noise and I am applying a window function and averaging the FFT^2? This is the code I am working with on matlab:

    n = 1;
freq1 = 1000.475;
freq2 = 1020.475;
x = 1:0.0001:2pi;
signal = 0.5sin(freq1(2pi)x+pi/5)+0.3sin(freq2(2pi)x+pi/5); % + 0.2sin(1000(2pi)x);% + 0.5sin(1000.8(2pi)*x);
sizex = size(x);
windowSize = sizex(2);
hanningWindow = hanning(windowSize);
noise_amplitude = 10;
totalIterations = 1000;  % total number of iterations in the loop
percentComplete = 0.1;    % the percentage at which you want to display progress
% Calculate the increments for progress display
progressIncrements = percentComplete:percentComplete:1;
% Initialize sum variables for Y and Ypad
sumY = zeros(1, windowSize);
sumYpad = zeros(1, 5 * windowSize);
for i=0:totalIterations
% Check if the current iteration corresponds to a progress increment
if any(abs(i / totalIterations - progressIncrements) < eps)
    % Display progress to the command window
    fprintf('%.0f%% completed\n', i / totalIterations * 100);
end

% signal plus noise
noise = noise_amplitude * randn(size(x));
y = signal + noise;

ywindowed = y.* hanningWindow';

% Zero-padding by a factor of 4
ywindowed_padded = [ywindowed, zeros(1, 4*length(ywindowed))];

if i==0
    figure(n); n = n+1;
    plot(x, ywindowed);
    title('Time Domain - Windowed Signal');
end

Y = abs(fft(ywindowed)).^2;

% zero padding

Ypad = fft(ywindowed_padded).^2;

% Accumulate Y and Ypad values
sumY = sumY + Y;
sumYpad = sumYpad + Ypad;

end
fprintf("Calculation done\n");
%         two sides    average            bin width         PSD    window function
sumYpad = 2sumYpad / (totalIterations  (length(ywindowed))^2);%/0.25;
sumY = 2sumY / (totalIterations  (length(ywindowed))^2);%/0.25;
%sumY = sumY -(sumY(10000)-sumYpad(10000)); % get padded and not padded to same level
figure(n); n = n+1;
clf;
frequencies = linspace(0, 1/(2*(x(2)-x(1))), length(sumYpad)/2);
plot(frequencies, abs(sumYpad(1:length(sumYpad)/2)));
xlim([999; 1001.5]);
%title('Frequency Domain - FFT of Windowed Signal with Zero Padding (Factor of 4)');
xlabel('Frequency (Hz)');
title("Frequency spectrum - interpolation - zero padding (factor of 4)")
hold on;
frequencies2 = linspace(0, 1/(2*(x(2)-x(1))), length(Y)/2);
% interpolate
xx = linspace(0, 1/(2*(x(2)-x(1))), length(sumY)/2);
yy = sumY(1:length(sumY)/2);
factor = 4;
xi = (linspace(min(xx), max(xx), length(xx)*factor))';
yi_spline = (interp1(xx, yy, xi, 'spline'));
plot(frequencies2, sumY(1:length(sumY)/2));
xlim([999; 1001.5]);
%title('Frequency Domain - FFT of Windowed Signal');
plot(xi, yi_spline, 'ro', 'MarkerSize', 5);
xlim([999; 1001.5]);
hold off
legend("zero padded", "original", "interpolation")

Answer 1

A lot of questions you have can probably be answered by this recent conversation I had here. Generally speaking, the power spectral density (PSD) doesn't contain information about the power in the signal. The magnitude squared of the output of the FFT is an estimate of the PSD, not the power spectrum. This can be proved as follows:

\begin{equation} S(\omega) = \lvert X(\omega)\rvert^{2} = X(\omega)X^{*}(\omega) = \mathcal{F}\{x(t)*x^{*}(-t)\}=\mathcal{F}\{r(t)\}\end{equation}

If your signal is real, meaning you have a mirrored spectrum, then it is possible to recover an estimate of the bandlimited power by computing

\begin{equation}\frac{1}{\pi}\int_{\omega_{1}}^{\omega_{2}}S(\omega)d\omega\end{equation}

but whether or not this is exact I'm not 100% sure.

For normalizing your spectrums, it is easiest to do normalize everything to 0 dB. This can be done by 10*log10(Y/max(abs(Y))). If you are wanting to normalize everything so that the minimum value on the magnitude spectrum is 0, as is the case in the other plot, you can do Y-min(abs(Y)). If you're looking for a different normalization, please clarify and I'll do my best to try and figure out how you want it normalized.

Windowing will affect the scaling of the spectrum depending on the summation of the window. If you add the line of code hanningWindow = hanningWindow/sum(hanningWindow) after defining the window, this should normalize the scaling effect of the windowing without altering how the window affects the signal's frequency content.

Lastly, likely the reason why the bin for your non-zero padded signal is not on the correct frequency location is because, based on the number of FFT points you are using, a bin doesn't exist at that location. This is the purpose of zero-padding, to interpolate values to extra bins based on where the integer multiples of the fundamental frequency exist. The peak lies between two integer multiples of the fundamental frequency, so it will get split into the two nearest bins. For whatever reason, it is favoring the one bin heavily over the other for accumulating amplitude.

Hopefully this helps, let me know if you have any other questions!