overlap-save方法 matlab实现

http://blog.sina.com.cn/s/blog_7ae6c3a00100sfa1.html

信号长度16384,包含频率100Hz和250Hz两个分量;
滤波器128点FIR低通,100Hz通过,250Hz被过滤.

load h; % 128-tap fir low pass filter: coefficients

N = 16*1024;
f1 = 100; % frequency of signal 1
f2 = 250; % frequency of signal 2
Fs = 2000; % sample rate

x1 = sin(2pif1*(0:N-1)/Fs); % signal 1
x2 = sin(2pif2*(0:N-1)/Fs); % signal 2

x = 0.8x1 + 0.7x2; % signal constructed by sig 1 and 2
%{
plot(x(1:256));
X1 = fft(x(1:256));
mag_X1 = abs(X1).^2
plot(mag_X1);
%}

M = length(h); % M = 128
N1 = 1024; % length of sub-segment
Lx = length(x); % Lx = 16384
L = N1-M; % L = 1024 - 128 = 896
hn = [h, zeros(1,N1-M)]; % padding zeros to h to get N1 length
POST_PAD = N1 - ((M+Lx) - floor((M + Lx)/L)L);
xn = [zeros(1,M), x, zeros(1,POST_PAD)];% zero padding at begin and end
K = floor((Lx+M)/L)+1; % K = 19
Y = zeros(K,N1); % Y = matrix 19x1024
for k=0:K-1
xk = xn(kL+1:k*L+N1); % each segment with M overlap
y = circconv(xk,hn,N1); % circle convolution
Y(k+1,:) = y; % a row of Y
end

Y = Y(:,M+1:N1)’; % discard first M points of each row then ’
y = (Y( : ))’; % convert to each column to a whole row.
y = y(1:Lx); % discard last POST_ADD points

plot(y(1:256)); % only display signal 1, signal 2 is filtered.
Y1 = fft(y(1:256)); % frequency domain
mag_Y1 = abs(Y1).^2; % magnitude
plot(mag_Y1); % display magnitude-frequency

//

用到的子程序1:

function [y,n] = circconv(x1,x2,N)
% implement circle convolution of x1 and x2
% x1 and x2 are input sequences of n1, n2 <= N
% y is output sequence
% N is size of circle buffer
% mehod: y(n) = sum(x1(m)*x2((n-m)mod N))
% check for length of x1

if (length(x1)>N |length(x2)>N)
error(‘length of input sequences is over range.’);
end
x1 = [x1 zeros(1,N-length(x1))];
x2 = [x2 zeros(1,N-length(x2))];
m = 0:N-1;
x2 =x2(mod(-m,N)+1);
H = zeros(N,N);
for n=1:1:N
H(n,:) = circshift_test(x2,n-1,N);
end

y = x1*H’;

用到的子程序2:

function y = circshift_test(x,m,N)
% circle shift of sequences
% Usage y = circshift_test(x,m,N)
% x: sequence to be circle shifted
% m: number to be shift to right
% N: length of circle shift
% y: output sequences

if (length(x)>N)
error(‘input length is over range.’);
end
x = [x zeros(1,N-length(x))];
n = 0:N-1;
n = mod(n-m,N);
y = x(n+1);