Hanning Window Energy Density Value

Good afternoon (code below figure),

I am processing some velocity data for a turbulent flow and want to determine the energy spectra to determine the shedding frequency (the picture below is from a test section that does not have a peak at 1.35Hz). I want to get rid off the noise. Format of the input data is:

time = 0, 0.05, 0.1, 0.15.

data = 1,1.2,1.04,0.95, ...

My colleague created a code that looks like this (fs = 200) and results in the "noisy" data plot at about 10^-2. Here data input has 4 coloms for 4 different velocity measurements.

len = length(data);

psd = (abs(fft(data(:,:),len)).^2)/(len*fs); % Fourier transform (2nd argument = vector length, such that the DFT returns n-points). Amplitude squared for PSD estimate, as this equals power

psd = 2.*psd(2:len/2+1,:); % single-sided spectrum. (1:N/2+1) takes only one half of the spectrum; the other half is its reflection

f = fs/2*linspace(0,1,len/2+1).'; % Nyquist frequency of the signal = fs/2

f = f(2:end); % remove zero Hz component

When I try and create the same with a hanning window; the magnitude is way off. Even multiplying it with 2 for amplitude correction does not help. Am I not seeing something?

Or how can I add the window to the code above? Is there a simple way?

Current code:

figure(1)

ylabel('Power Density Spectrum')

xlabel('Frequency (Hz)');

u_velocity=u_i(:,1); %Colum 1 is for u-velocity

mean_u = mean(u_velocity);

data = u_velocity-mean_u; %Velocity fluctuations around mean

chunk_size = 5000; % Size of each processing chunk

num_chunks = length(data) / chunk_size;

psd_combined = zeros(1, chunk_size / 2);

for i = 1:num_chunks

chunk_start = (i - 1) * chunk_size + 1;

chunk_end = i * chunk_size;

% Extract the current chunk of data

chunk = data(chunk_start:chunk_end);

% Apply the Hanning window

window = hanning(length(chunk));

windowed_data = chunk .* window';

% Perform FFT

fft_result = fft(windowed_data);

% Calculate Power Spectral Density (PSD) for this chunk

psd_chunk = abs(fft_result).^2;

% Combine the PSD results

psd_combined = psd_combined + psd_chunk(1:length(psd_chunk)/2);

end

psd_combined = psd_combined*10;%Correct Window

% Frequency axis

sample_rate = 200; % Adjust the sample rate accordingly

frequencies = linspace(0, sample_rate/2, length(psd_combined));

% Plot the PSD

semilogy(frequencies, psd_combined);

grid on;

Answer 1

[William Rose]

@Selina,

The normalization of the power spectrum is tricky. For example,there is the power spectrum and the power spectral density. Applyng a window reduces power somewhat. When you divide the signal into chunks and verage the spectra from the different chunks, this can introduce other normalization issues. Therefore I would not worry too much about the normalization. JUst compute your spectra in a consistent way, and then compare spectra (computed the same way) under different experimental conditions, or at different places in your system, to learn whatever it is you are trying to learn.

Your approach of dividng the signal into chunks and averaging the spectra from the different chunks is definitely a good one. And so is applying a window to each chuunk. This will reduce the noisiness of your spectra (i.e. the variance of your spectral esitmates) considerably.

I recommend using

pxx = pwelch(x,window,noverlap,nfft)

and I recommend using half-overlapped windows. This is because windowing tapers the signal to zero at the end of each chunk, so you are kind of throwing away data whn you do this. By using half-overlapped windows, you capture the parts that would otherwise be lost.

[William Rose]

@Selina,

[edit: add ";" on 'colors=...' line]

Since you did not provide data, I will use a random Gaussian signal, sampled at 200 Hz like your, and with duration 300 seconds, which is about the same as yours, I think, based on the lowest frequency in the spectrum with no windowing. It also appears from your spectra that you used chunks with duration 25, 10, 5, 2.5, and 1.25 seconds. In the plot below, note that the mean spectrum power is independent of the chunk size chosen - even when the chunk size = 300 s = full length of the signal.

fs=200;                 % sampling rate (Hz)
Ttot=300;               % duration (s)
Ns=fs*Ttot;             % number of samples
Tchunk=[300, 25, 10, 5, 2.5, 1.25];  % chunk duration (s)
Ncs=length(Tchunk);     % number of chunk sizes
Nchunk=Tchunk*fs;       % samples per chunk
x=randn(1,Ns);          % random Gaussian signal
colors=['k','r','g','b','c','m'];
for i=1:Ncs
    window=hann(Nchunk(i));
    noverlap=Nchunk(i)/2;  % for half-overlapped chunks
    [pxx,f] = pwelch(x,window,noverlap,Nchunk(i),fs); % compute spectral estimate
    % plot spectral estimate
    loglog(f,pxx,'-','Color',colors(i))
    hold on
end
grid on; xlabel('Frequency (Hz)'); ylabel('Power')
legend('300','25','10','5','2.5','1.25','location','southwest');

Since this is white noise, the expected spectrum will be perfectly flat.

The results above demonstrate how the use of many short chunks reduces the spectral variance. However, the lowest frequency in the spectrum is higher when you use short chunks, and the frequency resolution of the spectra is less good when you use short chunks.

[William Rose]

@Selina,

As long as you use a consistent method to compute the spectra, you should be able to identify the shedding frequency from the shape of the spectra, and the normalization will not affect your ability to do this.

The code below shows how to normalize the FFT of x, and adjust for the effect of the window, in order to get a power spectral density that is identical to the PSD returned by pwelch(), using one big chunk. I don't recommend using one big chunk. Use a lot of shorter chunks, to get a smoother spectrum. This just shows how the normalization of the PSD works.

fs=200; Ttot=300; Ns=fs*Ttot; % sampling freq (Hz); duration (s); number of samples
df=1/Ttot;                  % spectral resolution (Hz)
x=randn(1,Ns);              % random Gaussian signal
xwin=x.*hann(Ns)';          % windowed signal
sf=mean(hann(Ns).^2);       % power scale factor for Hann window
X2=abs(fft(xwin)/Ns).^2;    % raw power spectrum
X2=[X2(1),2*X2(2:Ns/2+1)];  % singled-sided power spectrum
psd=X2/df;                  % convert power spectrum to power spectral density
psd=psd/sf;                 % adjust for power loss due to window
f=[0:Ns/2]*df;              % frequencies for single-sided spectrum
% compute PSD with pwelch, using Hann window and one big chunk
[pxx,fw] = pwelch(x,hann(Ns),Ns/2,Ns,fs); % compute spectral estimate
% plot results
subplot(211), loglog(f,psd,'-r'); title('P.S.D. using FFT'); grid on
subplot(212), loglog(fw,pxx,'-r'); title('P.S.D. using pwelch'); grid on

The two spectra above are identical.

[Selina]

Thank you for your response! I tried running it with my data (attached as u_i but I am only interested in the first colum) and the two methods described in your comment showed the same result - here you can see the peak I was speaking about at about 1.3Hz. I see that you have been utilizing the built in functions. In order to adjust the window length, I tried changing NS to 10,000. However, it only resulted in error messages as the length x was no longer equal to NS (line xwin). This is why I used the loop beforehand. Is it possible to achieve this for both methods somehow but simpler then I tried to to do this? I am looking to clear up the signal but make sure the peak can be seen.

Thank you again! This was very helpful already :)

Edit: I assume the Ns/2 referred to half overlap?

[William Rose]

@Selina,

load('u_i.mat');
x=detrend(u_i(:,1));        % x=column 1, detrended
fs=200;                     % sampling freq (Hz)
Ns=length(x);               % number of samples
Tchunk=10;                  % chunk duration (s)
Nchunk=Tchunk*fs;           % samples per chunk
% compute PSD with pwelch, using Hann window
[pxx,f] = pwelch(x,hann(Nchunk),Nchunk/2,Nchunk,fs); % compute spectral estimate
% plot results
loglog(f,pxx,'-r'); title('Power Spectral Density'); grid on
xlabel('Frequency (Hz)'); ylabel('P.S.D. (power/Hz)');

I chose Tchunk=10 s, to obtain a good tradeoff between smoothing the spectrum while still retaining sufficient frequency resolution to resolve the peak at about 1.3 Hz. You can change Tchunk to see if you like the spectrum better with some other value, such as Tchunk=2.5, 5, 25, 50, etc.

I removed the linear trend from x(t) before computing the spectral estimate. This done to remove power at f=0, which is annoying on the plot, and which can leak into the nearby non-zero low frequencies, if it is not removed.

Yes, I chose noverlap=Nchunk/2 in order to obtain half-overlapped chunks.