graphics_toolkit gnuplot
pkg load signal
darkgreen = [33 150 33]/256;
M=7200; % big number, divisible by 8 and 9
% Generate M+1 samples of a triangle window
window = triang(M+1)'; % row vector
N=8; % actual window size, in "hops"
% Sample the triangle.
% Scale the abscissa. 0:M samples --> 0:9 "hops", and take 9 symmetrical hops, from .5 to 8.5
sam_per_hop_9 = M/9;
symmetric9 = window(1+(.5:8.5)*sam_per_hop_9);
periodic8 = symmetric9(1:8);
periodic_summation = [symmetric9(1)+symmetric9(N+1) symmetric9(2:N)];
% Re-scale the abscissa. 0:M samples --> 0:8 "hops", and take 8 symmetrical hops, from .5 to 7.5
sam_per_hop_8 = M/8;
symmetric8 = window(1+(.5:7.5)*sam_per_hop_8);
%------------------------------------------------------------------
% Compare windows based on their processing gain (PG) (Harris,1978,p 56,eq 15), because the ENBW
% formula allows values less than one "bin" (for some windows) when used with a 9-point periodic
% summation and an 8-point DFT. That actually makes sense, because a bandwidth of 1.1 (for instance)
% measured in 1/9-width bins is only 0.98 measured in 1/8-width bins. But values less than one
% are not customary, which could cause distrust.
PG_symmetric8 = sum(symmetric8)^2/sum(symmetric8.^2) % 6.0968
PG_periodic8 = sum(periodic8)^2 /sum(periodic8.^2) % 6.4529
PG_symmetric9 = sum(symmetric9)^2/sum(symmetric9.^2) % 6.7529
% Also note that the correct incoherent "power" formula for the
% periodic_summation window is sum(symmetric9.^2),
% not sum(periodic_summation.^2), because
% E{(h(1)·X(1) + h(9)·X(9))^2} = (h(1)^2 + h(9)^2)·E{X^2},
% not (h(1)^2 + 2·h(1)·h(9) + h(9)^2)·E{X^2}.
%------------------------------------------------------------------
% Plot the points
figure("position", [1 1 700 400])
plot(0:7, symmetric8, "color","red", ".", "markersize",14)
hold on
plot(8, symmetric9(9), "color","green", ".", "markersize",14)
plot(0:7, periodic8, "color","blue", ".", "markersize",14)
% Connect the dots
hops = (0:M)/sam_per_hop_9 -.5;
plot(hops, window, "color","blue") % periodic
hops = (0:M)/sam_per_hop_8 -.5;
plot(hops, window, "color","red") % symmetric
xlim([-.5 8.5])
set(gca, "xgrid","on")
set(gca, "ygrid","on")
set(gca, "ytick",[0:.25:1])
set(gca, "xtick",[0:8])
text(3.98, 0.69, 'L = 8 \rightarrow', "color","red", "fontsize",12)
text(5.25, 0.74, '\leftarrow L = 9', "color","blue", "fontsize",12)
title("Triangular window functions", "fontsize",14, "fontweight","normal")
xlabel('\leftarrow n \rightarrow', "fontsize",14)
% After this call, the cursor units change to a normalized ([0,1]) coordinate system
annotation("textarrow", [.76 .911], [.177 .2], "color",darkgreen,...
"string",{"discarded OR added to value at n=0 ";...
" (periodic summation)"}, "fontsize",10,...
"linewidth",1.5, "headstyle","vback1", "headlength",5, "headwidth",5)
%Now compute and plot the DTFTs and DFTs
M=64*N; % DTFT size
dr = 80; % dynamic range (decibels)
%------------------------------------------------------------------
% DTFT of symmetric window
H = abs(fft([symmetric8 zeros(1,M-N)]));
H = fftshift(H);
H = H/max(H);
H = 20*log10(H);
H = max(-dr,H);
x = N*[-M/2:M/2-1]/M;
figure("position", [1 1 700 400])
plot(x, H, "color","red", "linewidth",1);
hold on
ylim([-dr 0])
% Compute a DFT to sample the DTFT 8 times
H = abs(fft(symmetric8));
H = fftshift(H);
H = H/max(H);
H = 20*log10(H);
H = max(-dr,H);
plot(-N/2:(N/2-1), H, "color","red", ".", "markersize",14)
%------------------------------------------------------------------
% DTFT of periodic window
H = abs(fft([periodic8 zeros(1,M-N)]));
H = fftshift(H);
H = H/max(H);
H = 20*log10(H);
H = max(-dr,H);
plot(x, H, "color","blue", "linewidth",1);
% Compute a DFT to sample the DTFT 8 times
H = abs(real(fft(periodic8))); % real() is redundant... just to illustrate a point
H = fftshift(H);
H = H/max(H);
H = 20*log10(H);
H = max(-dr,H);
plot(-N/2:(N/2-1), H, "color","blue", ".", "markersize",14)
%------------------------------------------------------------------
% DTFT of a 9-sample symmetric window
H = abs(fft([symmetric9 zeros(1,M-N-1)]));
H = fftshift(H);
H = H/max(H);
H = 20*log10(H);
H = max(-dr,H);
plot(x, H, "color","green", "linewidth",1);
% Compute a DFT to sample the DTFT only 8 times.
H = abs(real(fft(periodic_summation))); % real() is redundant... just to illustrate a point
H = fftshift(H);
H = H/max(H);
H = 20*log10(H);
H = max(-dr,H);
plot(-N/2:(N/2-1), H, "color","green", ".", "markersize",14)
set(gca,"XTick", -N/2:N/2-1)
grid on
text(1.41, -18.8, {'\leftarrow DTFT';" symmetric 8"}, "color","red",...
"fontsize",10, "fontweight","bold")
set(gca,"XTick", -N/2:N/2-1)
grid on
ylabel("decibels", "fontsize",14)
xlabel("DFT bins", "fontsize",12, "fontweight","bold")
title('"Spectral leakage" from three triangular windows', "fontsize",14, "fontweight","normal")
% After this call, the cursor units change to a normalized ([0,1]) coordinate system
annotation("textarrow", [.132 .132], [.74 .6],...
"color", "blue", "string", {" DTFT";"periodic8"}, "fontsize",10,...
"linewidth",1, "headstyle","vback1", "headlength",5, "headwidth",5)
annotation("textarrow", [.28 .28], [.74 .613],...
"color", darkgreen, "string", {" DTFT";"symmetric9"}, "fontsize",10,...
"linewidth",1, "headstyle","vback1", "headlength",5, "headwidth",5)
annotation("arrow", [.524 .417], [.565 .752],...
"color", darkgreen,...
"linewidth",1, "headstyle","vback1", "headlength",5, "headwidth",5)
annotation("arrow", [.524 .632], [.565 .752],...
"color", darkgreen,...
"linewidth",1, "headstyle","vback1", "headlength",5, "headwidth",5)
annotation("textarrow", [.524 .524], [.525 .565], "color",darkgreen, "fontsize",10,...
"string",{" DFT 8";"periodic summation"},...
"linewidth",1, "headstyle","ellipse", "headlength",3, "headwidth",3)
% annotation("arrow", [.524 .311], [.565 .565], "linestyle", "--",...
% "color", darkgreen,...
% "linewidth",1, "headstyle","vback1", "headlength",5, "headwidth",5)
% annotation("arrow", [.524 .738], [.565 .565], "linestyle", "--",...
% "color", darkgreen,...
% "linewidth",1, "headstyle","vback1", "headlength",5, "headwidth",5)