Matlab for Spectrum Analysis Windows

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;

Matlab listing: dpssw.m

function [w,A,V] = dpssw(M,Wc);

% DPSSW - Compute Digital Prolate Spheroidal Sequence window of
% 	  length M, having cut-off frequency Wc in (0,pi).

k = (1:M-1);
s = sin(Wc*k)./ k;
c0 = [Wc,s];
A = toeplitz(c0);
[V,evals] = eig(A); % Only need the principal eigenvector
[emax,imax] = max(abs(diag(evals)));
w = V(:,imax);
w = w / max(w);

Sign in