フィルタ設計

過去の記事で,ESP32-S3でSDRを構築するためのコードを示してきたが,その中でフィルタの係数については具体的なものは示していなかったので,フィルタ係数を求める方法についてまとめておく.

ツールとしては GNU Octave を用いる.

以下に示すスクリプトを GNU Octave で実行すればフィルタ係数が記述されたヘッダファイルが生成される.

 

1.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]);

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト /  変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト /  変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト /  変更 )

%s と連携中