Multirate digital signal processing

% % m-file to illustrate simple interpolation and
% decimation operations (Program 9B.1, p641)
% File name: prog9b1.m
% An Illustration of interpolation by a factor of 4 %
Fs=1000; % sampling frequency
A=1.5; % relative amplitudes
B=1;
f1=50; % signal frequencies
f2=100;
t=0:1/Fs:1; % time vector
x=A*cos(2*pi*f1*t)+B*cos(2*pi*f2*t); % generate signal
y=interp(x,4); % interpolate signal by 4
stem(x(1:25)) % plot original signal
xlabel('Discrete time, nT ')
ylabel('Input signal level')
figure
stem(y(1:100)) % plot interpolated signal.
xlabel('Discrete time, 4 x nT')
ylabel('Output signal level')
==========================================

% % m-file to illustrate simple interpolation and
% decimation operations (Program 9B.2, p644).
% File name: prog9b2.m
% An Illustration of sampling rate changes using upfirdn by a factor of 4
%
Fs=1000; % sampling frequency
A=1.5; % relative amplitudes
B=1;
f1=50; % signal frequencies
f2=100;
t=0:1/Fs:1; % time vector
x=A*cos(2*pi*f1*t)+B*cos(2*pi*f2*t); % generate signal
y=resample(x,4,1); % interpolate signal by 4
stem(x(1:25)) % plot original signal
xlabel('Discrete time, nT ')
ylabel('Input signal level')
figure
stem(y(1:100)) % plot interpolated signal.
xlabel('Discrete time, 4 x nT')
ylabel('Interpolated output signal level')
y1=resample(y,1,4);
figure
stem(y1(1:25)) % plot decimated signal.
xlabel('Discrete time, nT')
ylabel('Decimated output signal level')

============================================
function moptimum
%Program moptimum is for designing I-stage optimum decimator
%or interpolator (I=1,2,3 or 4). The program computes the decimation
%factors, filter characteristics, and decimator efficiencies

%The following parameters must be provided by the user:
% Fs - input sampling frequency
% M - overall decimation factor
% fp - passband edge frequency
% dp - overall passband deviation in +ve dB
% ds - overall stopband deviation in +ve dB

clear all;
Fs = 96000; % sampling frequency in Hz
fp = 450; % passband edge frequency in Hz
dp = 0.0864; % overall passband deviation in +ve dB
ds = 60; % overall stopband deviation in +ve dB
M = 96; % overall decimation factor

EvalNStageDecimator(Fs,fp,dp,ds,M); % evaluate single stage decimator
for i=2:4 % evaluate all possible 2-, 3- and 4-stage decimators
R = GetFactors(M,i);
for j=1:size(R,1);
EvalNStageDecimator(Fs,fp,dp,ds,R(j,:));
end
end

===============================================
% % m-file for working out the computational
% complexities of a 2-stage decimator % Program name: mrate-ex1.m %
dp=0.01;
ds=0.01;
dp1=dp/2;
ds1=ds;
fp=5000;
Fi=1536000; F0=12000;
M=Fi/F0; M1=16; M2=8;
F1=Fi/M1; F2=F1/M2;
fs1=F1-(Fi/(2*M)); fs2=F2-(Fi/(2*M));
df1=(fs1-fp)/Fi; df2=(fs2-fp)/F1;
NUM= -10*log10(dp1*ds1)-13;
N1=(NUM/(14.6*df1)); N2=(NUM/(14.6*df2));
MPS=(N1*F1 + N2*F2); TSR = N1+N2;
M1
M2
fs1
fs2
df1
df2
N1
N2
MPS
======================================
function DecimationFactors = GetFactors(M,n)
% The function GetFactors finds all possible decimation factors for
% 2-, 3-, and 4-stage decimation. M is the
% overall decimation factor and n is the number of stages

p = floor(M/(2^(n-1)));
m = 1;
for i=2:p
for j=2:p
if n==2&i*j==M
R(m,1) = i; % get the 2-stage decimator factors
R(m,2) = j;
m = m + 1;
elseif n>2
for k=2:p
if n==3&i*j*k==M
R(m,1) = i; % get the 3-stage
R(m,2) = j; % decimator factors
R(m,3) = k;
m = m + 1;
elseif n>3
for l=2:p
if i*j*k*l==M
R(m,1) = i; % get the 4-stage
R(m,2) = j; % decimator factors
R(m,3) = k;
R(m,4) = l;
m = m + 1;
end
end
end
end
end
end
end
R = fliplr(sort(R')'); % sort the decimation factor vectors
z = zeros(1,size(R,2));
k = 1;
for i=1:size(R,1) % reject the redundancies
for j=i+1:size(R,1)
if R(i,:)==R(j,:)
R(j,:) = z;
end
end
if R(i,:)~=z
DecimationFactors(k,:) = R(i,:);
k = k + 1;
end
end

==========================================================
function decimation
%program performs multi-stages of decimation on data in a user-specified data file
%enough data must be guaranteed when using a large overall decimation fator

clear all;

r = [2 2 3 2]; % decimation factor array for different stages
FIR = 0; % 1 - use FIR filter, 0 - use IIR filter
n = 0; % order of IIR filter or FIR filter length
% n=0 for 30 points FIR filter or 8th order Chebyshev type I LPF filter
in = fopen('decimation.dat','r') ;
[x,count] = fscanf(in,'%g',inf); % data to be decimated
y = x;
fclose(in);
for i=1:length(r)
if n==0 %decimating data use default filter
if FIR
y = decimate(y,r(i),'fir');
else
y = decimate(y,r(i));
end
else
if FIR
y = decimate(y,r(i),n,'fir');
else
y = decimate(y,r(i),n);
end
end
end
plot(1:count,x,1:prod(r):count,y(1:length(y)),'o');

===========================================================
function EvalNStageDecimator(Fs,fp,dp,ds,R)
format long;
a1 = 0.005309;
a2 = 0.07114;
a3 = -0.4761;
a4 = -0.00266;
a5 = -0.5941;
a6 = -0.4278;
a7 = 11.01217;
a8 = 0.51244;
dp = 10^(dp/20.0)-1; ds = 10^(-ds/20.0);
Ftemp = Fs;

dp = dp/length(R);
MPS = 0; TSR = 0;
fprintf('stage\tfactor\tFi\tfp\tfs\tdp\tds\tN\n');
for i=1:length(R) % n=size(R,1) possible k-stages decimators
F = Ftemp/R(i);
fs = F - Fs/2/prod(R);
df = (fs - fp)/Ftemp;
Ftemp = F;
N = round((log10(ds)*(a1*log10(dp)^2+a2*log10(dp)+a3)+...
a4*log10(dp)^2+a5*log10(dp)+a6)/df-df*(a7+a8*(log10(dp)-log10(ds)))+1);
fprintf('%d\t%d\t%d\t%d\t%d\t%-7.4f\t%-7.4f\t%d\n',i,R(i),F,fp,fs,dp,ds,N);
MPS = MPS + N*F;
TSR = TSR + N;
end
fprintf('MPS = %d, TSR = %d\n\n',MPS,TSR);
format;

========================
function interpolation
%program performs multi-stages interpolation on data in a user-specified data file

clear all;

r = [2 1 1]; % decimation factor array for different stages
l = 4; % filter length
alpha = 0.5; % cut-off frequency
in = fopen('decimation.dat','r') ;
[x,count] = fscanf(in,'%g',inf); % data to be decimated
y = x;
fclose(in);
for i=1:length(r)
y = interp(y,r(i),l,alpha);
end
subplot(1,2,1);stem(x);subplot(1,2,2);stem(y);