Hanning Window Energy Density Value

Question

Selina 2023-11-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2047727-hanning-window-energy-density-value

评论： William Rose 2023-11-16

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;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

William Rose 2023-11-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2047727-hanning-window-energy-density-value#answer_1353152

@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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

William Rose 2023-11-15

编辑：William Rose 2023-11-15

在 MATLAB Online 中打开

@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 2023-11-15

编辑：William Rose 2023-11-15

在 MATLAB Online 中打开

@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 2023-11-15

编辑：Selina 2023-11-15

u_i.mat

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 2023-11-16

在 MATLAB Online 中打开

u_i.mat

@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.

请先登录，再进行评论。

Hanning Window Energy Density Value

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Hanning Window Energy Density Value

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论