《DSP using MATLAB》Problem 5.32

代码：

function [y] = ovrlpadd_v3(x, h, N)

%% Overlap-Add method of block convolution

%% ------------------------------------------------------

%% [y] = ovrlpadd(x, h, N)

%% y = output sequence

%% x = input sequence

%% h = impulse response

%% N = block length >= 2*length(h)-1

N = 2^(ceil(log10(N)/log10(2)))

Lx = length(x);                        % ML

L = length(h);                         % length of impulse response

Hk_FFT = fft(h, N);

M = floor((Lx)/(L))                   % number of bolck

Y = zeros(1, Lx+N-L)

% convolution with successive blocks

for m = 0:M-1

	xm = x(m*L+1 : m*L+L);                                     % length is L

	XMk_FFT = fft(xm, N);

	YMk_FFT = XMk_FFT .* Hk_FFT;

	ym = real(ifft(YMk_FFT))                                   % length is 2L-1

	%Y(m*L+1 : m*L+(2*L-1)) = Y(m*L+1 : m*L+(2*L-1)) + ym

	Y(m*L+1 : m*L+N) = Y(m*L+1 : m*L+N) + ym

end

y = Y;

　　

%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%%            Output Info about this m-file

fprintf('\n***********************************************************\n');

fprintf('        <DSP using MATLAB> Problem 5.32 \n\n');

banner();

%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

% -------------------------------------------------------------

%           Overlap-Add

% -------------------------------------------------------------

n1 = [0:8-1];

x1 = cos(pi*n1/500);

N1 = length(x1);

nh = [0:3];

 h = [1, -1, 1, -1];

 L = length(h);

M = N1/L;

N = 2*L-1;

y = ovrlpadd_v3(x1, h, N);            % FFT

　　运行结果：

原题中x(n)长4000，不好看图说明，这里假设长度只有8，便于上图说明道理。