Why is my hamming window and rectangular window values the same for sigma_omega when I make L = 21?

3 次查看(过去 30 天)
I have this code that produces the same value for both the hamming window and the rectangular window when I try to solve for sigma_omega for given L (L = 21, 51, and 101). For the sake of the code I will only be testing for L = 21.
Here is the equation for sigma_omega:
Here is the code I have:
% Set window length L = 21
L = 21;
% Generate rectangular and Hamming windows
rect_window = window(@rectwin, L);
hamming_window = window(@hamming, L);
% Calculate the frequency spread for each window
sigma_omega_rect = calc_freq_spread(rect_window);
sigma_omega_hamming = calc_freq_spread(hamming_window);
% Display the frequency spread for rectangular and Hamming windows for L = 21
fprintf('L = 21: sigma_omega (Rectangular) = %.4f\n', sigma_omega_rect);
L = 21: sigma_omega (Rectangular) = 58.0983
fprintf('L = 21: sigma_omega (Hamming) = %.4f\n\n', sigma_omega_hamming);
L = 21: sigma_omega (Hamming) = 58.0983
% Function to calculate frequency spread
function sigma_omega = calc_freq_spread(w)
L = length(w);
% Calculate DTFT (using FFT as an approximation)
N = 1024; % Length of the FFT to improve frequency resolution
W = fft(w, N); % Compute the FFT of the window
% Frequency axis (normalized frequency ω from -π to π)
omega_hat = linspace(-pi, pi, N);
% Magnitude squared of DTFT
W_mag_sq = abs(W).^2;
% Numerator: Sum of ω^2 * |W(e^jω)|^2 (approximating the integral)
numerator = sum((omega_hat.^2) .* W_mag_sq);
% Denominator: Sum of |W(e^jω)|^2 (approximating the integral)
denominator = sum(W_mag_sq);
% Calculate the frequency spread (σ_ω)
sigma_omega = sqrt(sum(numerator) / sum(denominator)); % Summing the whole thing to a scalar
end
My question is what is wrong with my code for producing the same output? The same exact output for both rectangular and hamming window is also the same 58.0983 when I change L to either 51 or 101. Please advise.
Thank you.

采纳的回答

Matt J
Matt J 2024-10-6
编辑:Matt J 2024-10-6
% Set window length L = 21
L = 21;
% Generate rectangular and Hamming windows
rect_window = window(@rectwin, L);
hamming_window = window(@hamming, L);
% Calculate the frequency spread for each window
sigma_omega_rect = calc_freq_spread(rect_window);
sigma_omega_hamming = calc_freq_spread(hamming_window);
% Display the frequency spread for rectangular and Hamming windows for L = 21
fprintf('L = 21: sigma_omega (Rectangular) = %.4f\n', sigma_omega_rect);
L = 21: sigma_omega (Rectangular) = 0.3635
fprintf('L = 21: sigma_omega (Hamming) = %.4f\n\n', sigma_omega_hamming);
L = 21: sigma_omega (Hamming) = 0.1685
% Function to calculate frequency spread
function sigma_omega = calc_freq_spread(w)
% Calculate DTFT (using FFT as an approximation)
N = 1024; % Length of the FFT to improve frequency resolution
w=w(:).'; %<---Matt J
omega_hat = ((0:N-1)-ceil((N-1)/2))*2/N*pi; %<---Matt J
W = fftshift(fft(w, N)); % <--- Matt J
% Magnitude squared of DTFT
W_mag_sq = abs(W).^2;
% Numerator: Sum of ω^2 * |W(e^jω)|^2 (approximating the integral)
numerator = sum((omega_hat.^2) .* W_mag_sq);
% Denominator: Sum of |W(e^jω)|^2 (approximating the integral)
denominator = sum(W_mag_sq);
% Calculate the frequency spread (σ_ω)
sigma_omega = sqrt(sum(numerator) / sum(denominator)); % Summing the whole thing to a scalar
end

更多回答(1 个)

Paul
Paul 2024-10-6
编辑:Paul 2024-10-6
Another approach is to integrate the DTFT directly. The DTFT can be evaluted with freqz
% Set window length L = 21
L = 21;
% Generate rectangular and Hamming windows
rect_window = window(@rectwin, L);
hamming_window = window(@hamming, L);
dtft = @(x,w) freqz(x,1,w);
sigma_w = @(x) sqrt(integral(@(w) w.^2.*abs(dtft(x,w)).^2,-pi,pi)/integral(@(w) abs(dtft(x,w)).^2,-pi,pi));
sigma_w(rect_window)
ans = 0.3635
sigma_w(hamming_window)
ans = 0.1685
Taking advantage of Parseval for the denominator:
sigma_w = @(x) sqrt(integral(@(w) w.^2.*abs(dtft(x,w)).^2,-pi,pi)/(2*pi*sum(abs(x).^2)));
sigma_w(rect_window)
ans = 0.3635
sigma_w(hamming_window)
ans = 0.1685

产品


版本

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by