Blackman Window Example

Below is the Matlab script for creating Figures 2.6 and 2.7 in §2.5.4.

% Illustrate zero-phase zero-padding around a Blackman window

% Analysis parameters:
M = 31;         % Window length
N = 64;         % FFT length (zero padding factor = N/M)
Mo2 = (M-1)/2;  % Shorthand
dBCut = -100;   % Clip dB at this level

% Signal parameters (synthetic sinusoid):
wxT = pi/exp(1);% Sinusoid frequency in rad/sample
A = 1;          % Sinusoid amplitude
phix = pi/3;    % Sinusoid phase

% Compute the signal x:
n = [0:N-1];    % time indices for sinusoid and FFT
x = A * cos(wxT*n+phix); % a sinusoid

% Compute a causal Blackman window:
% w = blackman(M); % Matlab signal processing toolbox, or
w = .42-.5*cos(2*pi*(0:M-1)/(M-1))+.08*cos(4*pi*(0:M-1)/(M-1));

% Apply the window:
xw = w .* x(1:M);
% Note that we have "skipped" the first Mo2 samples of x
% since the center of the window is its proper time origin.

% Form the zero-padded FFT input buffer.
% Note how the negative-time portion goes on the right:
xwzp = [xw(Mo2+1:M),zeros(1,N-M),xw(1:Mo2)];
wzp = [w(Mo2+1:M),zeros(1,N-M),w(1:Mo2)];

figure(1);
subplot(1,1,1);  % force a clear
subplot(2,1,1); 
n = -Mo2:Mo2;    % time axis for plot
stem(n,xw,'ok'); % windowed data
hold on;
plot(n,w,':k');  % window
plot(n,zeros(1,M),'-k');  % zero line
title('Blackman Windowed Sinusoid');
xlabel('Time (samples)'); 
ylabel('Amplitude');
text(-8,1,'(a)'); 
hold off;

subplot(2,1,2); 
n = 0:N-1;     % FFT buffer time axis 
stem(n,xwzp,'ok'); 
hold on;
plot(n,zeros(1,N),'-k');  % zero line
plot(n,wzp,':k');         % window
plot([N/2,N/2],[-1,1],'--k');  % dividing line
text(N/4,0.2,'positive time'); 
text(N/2+N/8,0.2,'negative time'); 
xlabel('Time (samples)'); 
ylabel('Amplitude');
text(-8,1,'(b)'); 

saveplot('../eps/zpblackmanT.eps');
hold off;

figure(2);
subplot(1,1,1); % force a clear

% Now show the window transform:
Xwzp = fft(xwzp);             % Blackman window transform
spec = 20*log10(abs(Xwzp));   % Spectral magnitude in dB
spec = spec - max(spec);      % Normalize to 0 db max
spec = max(spec,dBCut*ones(1,N)); % clip to dBCut dB

subplot(2,1,1); 
fni = [0:1.0/N:1-1.0/N];   % Normalized frequency axis
bins = [0:N-1];            % Bin axis
Xwzpa = abs(Xwzp)
stem(bins,Xwzpa,'ok'); 
hold on;
plot([N/2,N/2],YLim,'--k'); % dividing line
text(5,7,'positive frequencies'); 
text(45,7,'negative frequencies'); 
%zpblackmanplot3d
axis([0,N-1,YLim]); 
xlabel('Frequency (bins))'); 
ylabel('Magnitude (linear)');
text(-8,8,'(a)'); 

% Replot interpreting upper bin numbers as (-) frequencies:
subplot(2,1,2);
sbins = bins - N/2;
nh = N/2;
sXwzpa = [Xwzpa(nh+1:N),Xwzpa(1:nh)]; % see also fftshift()
%fninf = fni - 0.5;               % Plot (-) freqs on the left
%plot(fninf,sXwzpa,'-k'); axis([-0.5,0.5,dBCut,10]); grid;
%plot(fni,Xwzpa,'-k'); grid;
stem(sbins,abs(sXwzpa),'ok');
axis([-N/2,N/2-1,YLim]); 
hold on;
plot([0,0],YLim,'--k');  % dividing line
text(-20,7,'negative frequencies'); 
text(5,7,'positive frequencies'); 
xlabel('Frequency (bins))'); 
ylabel('Magnitude (linear)');
text(-48,8,'(b)'); 
saveplot('../eps/zpblackmanF.eps');
hold off;