clear all; close all; clc;
% filter specification:
% sampling-rate conversion factor M=5.
% passband edge at OMEGA_p = 0.09pi. stopband edge at OMEGA_s=pi/5.
% passband and stopbands ripples are given in decibels with ap=0.1dB,
% as=60dB.
% step 1: determining the input parameters for FIR filter desgin
ap = 0.1; as=60; Fp=0.09*pi/pi; Fs=pi/5/pi; % filter specifications
dev=[(10^(ap/20)-1)/(10^(ap/20)+1), 10^(-as/20)]; % passband and stopband ripples
F = [Fp,Fs]; % cutoff frequencies
A = [1,0]; % desired amplitudes
[Nord, Fo, Ao, W] = firpmord(F,A,dev); % estimating the filter order
% increase Nord to match stopband ripple requirement
Nord = Nord+4;
h = firpm(Nord, Fo, Ao, W); % filter coefficients
plot(h);
figure(1);
freqz(h,1,1024) % computes and plots the frequency response
axis([0 1 -80.5 1]);
figure(2);
subplot(2,1,1); impz(h,1) % plots the impulse response
subplot(2,1,2); zplane(h,1) % plots pole/zero
fs = 48e3; % sample rate
f1 = 1e3; % signal1: 1kHz
f2 = 16e3; % signal2: 16kHz
N = 8192; % signal length
f = (0:N/2-1)/N*fs;
x = 0.5*sin(2*pi*(0:N-1)*f1/fs)+0.5*sin(2*pi*(0:N-1)*f2/fs); % input signal
y = filter(h,1,x); % direct filtering using matlab function
block_size = 128; % split input signal into consequent 128-points block
filter_size = length(h); % fir filter tap
N_b = N/block_size; % how many blocks of total input signal
x_f = zeros(1,block_size+filter_size-1); % block_size + history_size
y_f = zeros(1,N_b*block_size); % signal filterred
for n=1:N_b % for each 128-points block
x_f(filter_size:end) = x((n-1)*block_size+1:n*block_size); % get a new block
for i=1:block_size % for each point in block
sum = 0;
for k=1:filter_size
sum = sum+x_f(k+i-1)*h(k); % h(k) is symmetrical,
% if h(k) is not symmetrical, need to invert h
end
y_f((n-1)*block_size+i) = sum; % output sample filterred
end
for i=1:filter_size-1
x_f(i) = x_f(i+block_size); % update history for each block
end
end
error = y-y_f; % error of different filtering procedure
figure(3);
plot(error); % error is close to zeros
从第28行开始看。
|
楼主
