過去の記事で，ESP32-S3でSDRを構築するためのコードを示してきたが，その中でフィルタの係数については具体的なものは示していなかったので，フィルタ係数を求める方法についてまとめておく．
ツールとしては GNU Octave を用いる．
以下に示すスクリプトを GNU Octave で実行すればフィルタ係数が記述されたヘッダファイルが生成される．
１．IIR LPF（8次エリプティック LPF）
% This file name: "IIR_LPF_design.m"
%-----------------------------------------------
% To generate "IIR_LPF.h",
% execute this script for MATLAB/GNU Octave.
%
% Jan. 5, 2023 by T. Uebo
%-----------------------------------------------
pkg load signal;
clear;
close all;
clc;
%-- Modify this section if you need ----
fs=96e3; %[Hz]
fc=5e3; %{Hz]
Rp=1; %[dB]
Rs=100; %[dB]
%----------------------------------------
% 8th LPF
ord=8;
% Elliptic LPF
[b,a]=ellip(ord, Rp, Rs, 2*fc/fs);
% Chebyshev LPF
%[b,a]=cheby1(ord, Rp, 2*fc/fs);
figure(1);
freqz(b,a,2048,fs);
[sos, g]=tf2sos(b,a);
g=g.^(1/(ord/2));
B0=g.*sos(1,1:3);
A0=sos(1,5:6);
B1=g.*sos(2,1:3);
A1=sos(2,5:6);
B2=g.*sos(3,1:3);
A2=sos(3,5:6);
B3=g.*sos(4,1:3);
A3=sos(4,5:6);
FID=fopen("IIR_LPF.h", "w");
fprintf(FID,"// %dth IIR Elliptic LPF\n" ,ord);
fprintf(FID,"// fs=%d[Hz]\n" ,fs);
fprintf(FID,"// fc=%d[Hz]\n" ,fc);
fprintf(FID,"// Ripple=%d[dB]\n", Rp);
fprintf(FID,"// Att=%d[dB]\n", Rs);
fprintf(FID,"// \n");
fprintf(FID,"// b0 + b1*Z^(-1) + b2*z^(-2) \n");
fprintf(FID,"// H(z)=--------------------------- \n");
fprintf(FID,"// 1 + a1*z^(-1) + a2*z^(-2) \n");
fprintf(FID,"// \n");
fprintf(FID,"// {b0, b1, b2, a1, a2}\n");
fprintf(FID,"float biquad0[] =\n");
fprintf(FID,"{%e,%e,%e,%e,%e};\n",[B0 A0]);
fprintf(FID,"float biquad1[] =\n");
fprintf(FID,"{%e,%e,%e,%e,%e};\n",[B1 A1]);
fprintf(FID,"float biquad2[] =\n");
fprintf(FID,"{%e,%e,%e,%e,%e};\n",[B2 A2]);
fprintf(FID,"float biquad3[] =\n");
fprintf(FID,"{%e,%e,%e,%e,%e};\n",[B3 A3]);
fclose(FID);
2. FIR LPF
% This file name: "FIR_LPF_design.m"
%-----------------------------------------------
% To generate "FIR_LPF.h",
% execute this script for MATLAB/GNU Octave.
%
% Jan. 5, 2023 by T. Uebo
%-----------------------------------------------
pkg load signal;
clear;
close all;
clc;
%-- Modify this section if you need ----
Nfir=256;
fs=16e3; %[Hz]
fc=3e3; %[Hz]
%----------------------------------------
b = fir1(Nfir-1,2*fc/fs, 'low',...
blackmanharris(Nfir) );
figure(1)
freqz(b,1,2048, fs)
figure(2)
plot(b);
FID=fopen("FIR_LPF.h", "w");
fprintf(FID, "float FIR_LPF[]={\n");
for k=[1:Nfir-1]
fprintf(FID,"%e,\n", b(k));
end;
fprintf(FID,"%e\n", b(Nfir));
fprintf(FID,"};\n");
fclose(FID);
3. 複素係数FIR BPF
% This file name: "CPLX_BPF_design.m"
%-------------------------------------------------
% To generate "CPLX_filter.h",
% execute this script for MATLAB/GNU Octave.
%
% Jan. 5, 2023 by T. Uebo
%-------------------------------------------------
pkg load signal
clear;
close all;
clc;
%--- Modify this section if you need --------
% Define Pass Band: f1(Hz) to f2(Hz)
f1=250;
f2=2850;
%Sampling frequency
fs=16e3;
%Number of Taps
Ndat=256;
%--Set a Window function---
%wf=ones(Ndat,1);
%wf=hanning(Ndat);
%wf=hamming(Ndat);
wf=blackman(Ndat);
%wf=blackmanharris(Ndat);
%wf=chebwin(Ndat);
%--------------------------------------------
%--- Calc. complex FIR coeff. ---------------
if (f1<f2)
fa=f1+fs/2;
fb=f2+fs/2;
else
fa=f2+fs/2;
fb=f1+fs/2;
end
ka=round( fa/(fs/Ndat) );
kb=round( fb/(fs/Ndat) );
fres=[ zeros(1,ka-1)...
ones(1,kb+1-ka)...
zeros(1,Ndat-kb)];
fres=[fres(Ndat/2+1:end) fres(1:Ndat/2)];
imp_res=ifft(fres);
Ich=[real(imp_res(Ndat/2+1:end))...
real(imp_res(1:Ndat/2))];
Qch=[imag(imp_res(Ndat/2+1:end))...
imag(imp_res(1:Ndat/2))];
Ich=Ich.*wf.';
Qch=Qch.*wf.';
%----Generate coeff. file ------------------
FID=fopen('CPLX_filter.h', 'w');
fprintf(FID, 'float CF_Re[%d]={\n' , Ndat);
for p=[1:Ndat-1]
fprintf(FID,'%e,\n', Ich(p));
end;
fprintf(FID,'%e\n', Ich(Ndat) );
fprintf(FID,'};\n');
fprintf(FID, 'float CF_Im[%d]={\n' , Ndat);
for p=[1:Ndat-1]
fprintf(FID,'%e,\n', Qch(p) );
end;
fprintf(FID,'%e\n', Qch(Ndat) );
fprintf(FID,'};\n');
fclose(FID);
%---- Display results ------------------------
figure(1); plot(Ich);
ylabel('FIR coeff. (Real Part)');
figure(2); plot(Qch);
ylabel('FIR coeff. (Imaginary Part)');
figure(3);
[hr,fr]=freqz(Ich+j*Qch,[1], Ndat*16,'whole',fs);
fr=fr-fs/2;
hr=[hr(Ndat*8+1:end); hr(1:Ndat*8)];
plot(fr, 20.*log10( abs(hr)) );
xlabel('Frequency(Hz)');
ylabel('Response(dB)');
grid on;
figure(4)
subplot(1,2,1)
plot(fr, 20.*log10( abs(hr)) );
xlabel('Frequency(Hz)');
ylabel('Response(dB)');
grid on;
axis([f1-500, f1+500, -80 10]);
subplot(1,2,2)
plot(fr, 20.*log10( abs(hr)) );
xlabel('Frequency(Hz)');
ylabel('Response(dB)');
grid on;
axis([f2-500, f2+500, -80 10]);
Tj Lab 