Hanning Window Energy Density Value

9 visualizzazioni (ultimi 30 giorni)
Selina
Selina il 15 Nov 2023
Commentato: William Rose il 16 Nov 2023
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;

Risposta accettata

William Rose
William Rose il 15 Nov 2023
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
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.
  4 Commenti
Selina
Selina il 15 Nov 2023
Modificato: Selina il 15 Nov 2023
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
William Rose il 16 Nov 2023
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.

Accedi per commentare.

Più risposte (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by