Verilogでアナログシュミレート？
状態方程式の一般形式は
です。制御系の技術者にはおなじみの式ですね。これを４次のルンゲクッタで計算させています。
呼び出し形式は
| $runge_kutta("ファイル名”、モード、Input_RealArray、Output_RealArray）； |
それでは、使い方のデモソースです。
Initialの$runge_kuttaで、遷移行列ファイルの指定、モードの指定、入出力のREALアレーを定義します。
次にAlways文で、毎クロック呼び出します。
`timescale 1ns / 1ns
//Apr.18.2005 rungekutta IF changed.
module runge_kutta_test;
parameter freq =900e3;
parameter dac_width=10;
real read_array[0:10];
real write_array[0:10];
real read_array2[14:24];
real write_array2[10:24];
real t,amp,ar;
reg clock=0;
reg [15:0] filt_out1,filt_out2;
reg [dac_width-1:0] dac_input;
real range;
reg Reset;
wire dac_out;
real ar;
always #1 clock=~clock;
always @(posedge clock) begin
t=$time()*1e-9;
amp=($sin(2*$M_PI()*t*freq)+1)*range;
dac_input=$rtoi(amp);
read_array[0]=dac_out*6283185.3;//cut off 1e6Hz degree2
//ar=$sin(2*$M_PI()*t*freq);
read_array2[14]=dac_out*6283185.3;//cur off 1e6Hz degree10
$runge_kutta("Ae6_2.mat",1);
$runge_kutta("Ae6.mat",1);
ar=dac_out*6283185.3;
// filt_out1=write_array[1]*range;//2degree 0+(2-1)=1
// filt_out2=write_array2[19]*range;//10degree 10+(10-1)=19
filt_out1=write_array[0]*range;//Apr.18.2005 2degree 0+(2-1)=1
filt_out2=write_array2[10]*range;//Apr.18.2005 10degree 10+(10-1)=19
//$display("filt_out2=%e %e",ar,write_array2[11]);
end
sigma_delta #(dac_width) dac1(dac_out, dac_input,clock, Reset);
initial begin
Reset=1;
range=((1 << dac_width)-2)/2;
$runge_kutta("Ae6_2.mat",0,read_array,write_array);
$runge_kutta("Ae6.mat",0,read_array2,write_array2);
#100 Reset=0;
#100000 $finish;
end
endmodule
module sigma_delta(out, DACin, Clk, Reset);
parameter dac_with=10;
input [dac_with-1:0] DACin; //
input Clk;
input Reset;
output out;
reg [dac_with-1+2:0] R1;
wire out=R1[dac_with-1+2];
wire [dac_with-1+2:0] delta= {R1[dac_with-1+2], R1[dac_with-1+2], {dac_with{1'b0}}};
wire [dac_with-1+2:0] adder1 = {2'b00,DACin} + delta;
wire [dac_with-1+2:0] adder2 = adder1 + R1;
always @(posedge Clk or posedge Reset)
begin
if(Reset) R1 <= {1'b1, 11'h000};
else R1 <=adder2;
end
endmodule
|
本例では、一次のシグマデルタDAコンバータでDA変換しています。シグマデルタコンバータの出力は1ビットのビットストリームになります。これを2次、および１０次のアナログバターワースLPFでアナログに戻してみます。
１０次フィルタの行列は次で与えられます。
ファイルとして読み込むのは、Aマトリクスだけです。残りは、Verilogソース内でプログラムします。ルンゲクッタの計算結果xはwrite_array2に入ります。write_array2は、LEFT_ADDRESS１０から始まり10次の行列ですので出力は19にはいることになります。
A =
1.0e+008 *
-1.2412 -0.6283 0 0 0 0 0 0 0 0
0.6283 0 0 0 0 0 0 0 0 0
0 0.6283 -1.1197 -0.6283 0 0 0 0 0 0
0 0 0.6283 0 0 0 0 0 0 0
0 0 0 0.6283 -0.8886 -0.6283 0 0 0 0
0 0 0 0 0.6283 0 0 0 0 0
0 0 0 0 0 0.6283 -0.5705 -0.6283 0 0
0 0 0 0 0 0 0.6283 0 0 0
0 0 0 0 0 0 0 0.6283 -0.1966 -0.6283
0 0 0 0 0 0 0 0 0.6283 0
B =
1.0e+007 *
6.2832
0
0
0
0
0
0
0
0
0
C =
0 0 0 0 0 0 0 0 0 1.0000
D =
0
下がシミュレーション波形です。DACの1ビット出力DAC_OUTがPPM変調されている様子がわかります。これをフィルタに通すと、filt_out1,filt_out2になります。2次と10次の違いはあまりわかりませんが、とりあえず、10次の方が次数が大きい分位相Delayが大きいことがわかります。
最後は、VPI ソースリストです。
// Mar.19.2004 www.sugawara-systems.com
// Linear System Simulator on Verilog2001 VPI
// Apr.18.2005 Bug fix address to only 1 output, change output port0
//
#include "runge_kutta.h"
#include "math.h"
# include "..\..\vpi\vpi\vpi_user.h"
#include <cassert>
#include <map>
#include <string>
extern int reset_lexor(const char * path);
extern bool parse_mat_data();
extern vector< vector<double*> *> *doubles_types;
//static runge_kutta *gloval_ptr=0;
static map <string , runge_kutta*> file_map;
runge_kutta::runge_kutta()
{
doubles_list=0;
continue_runge_kutta=false;
old_time=0;
zh=0;
mode=INITIAL;
}
runge_kutta::~runge_kutta()
{
}
bool runge_kutta::file_read(string file_name)
{
bool success=reset_lexor(file_name.c_str());
if(success){
success=parse_mat_data();
if (success) doubles_list=doubles_types;
else {
vpi_mcd_printf(1,"$runge_kutta: can not read coefficient file.\n");
}
}else {
vpi_mcd_printf(1,"$runge_kutta:can not open coefficient file.\n");
}
return success;
}
bool runge_kutta::check_and_dump_double_list()
{
vpi_mcd_printf(1,"coefficient file has read.\n");
rows=doubles_list->size();
assert(rows);
columns= (*doubles_list)[0]->size();
bool error=false;
if (rows ==columns){
vpi_mcd_printf(1,"Matrix is [%d][%d]\n",rows,columns);
for (unsigned i=0; i< rows;i++) {
if ((*doubles_list)[i]->size() !=columns) error=true;
vector <double> a;
for (unsigned m=0; m < (*doubles_list)[i]->size();m++){
double temp=*(*(*doubles_list)[i])[m];
a.push_back(temp);
vpi_mcd_printf(1,"A Matrix[%d][%d]=%f \n",i,m,temp);
}
A.push_back(a);
}
if (columns !=rows) error =true;
if (!error) {
//Initialize array
for (unsigned i=0; i<rows;i++){
zx.push_back(0.0);
zdx.push_back(0.0);
zq.push_back(0.0);
}
for (unsigned i=0; i< rows;i++) {
if (i==0) B.push_back(1.0);//Jan.5.2005
else B.push_back(0.0);//Jan.5.2005
if (i==rows-1) C.push_back(1.0);//Jan.5.2005
else C.push_back(0.0);//Jan.5.2005
}
D=0.0;//Jan.5.2005
}else {
vpi_mcd_printf(1,"$runge_kutta:can not read coefficient file.\n");
}
return !error;
}else {
// x=A nxn B nx1
// y=C 1xn D 1x1
// n+n+1+1
if (rows ==2*columns+2){
//vpi_mcd_printf(1,"Matrix is/*行列は*/[%d][%d]./*行列です。*/\n",rows,columns);
for (unsigned i=0; i< columns;i++) {
if ((*doubles_list)[i]->size() !=columns) error=true;
vector <double> a;
for (unsigned m=0; m < columns;m++){
double temp=*(*(*doubles_list)[i])[m];
a.push_back(temp);
vpi_mcd_printf(1,"A Matrix[%d][%d]=%f \n",i,m,temp);
}
A.push_back(a);
}
//B
for (unsigned i=0; i< columns;i++) {
double temp=*(*(*doubles_list)[i+columns])[0];
vpi_mcd_printf(1,"B Matrix[%d][0]=%f \n",i,temp);
B.push_back(temp);
}
//C
for (unsigned i=0; i< columns;i++) {
double temp=*(*(*doubles_list)[2*columns])[i];
vpi_mcd_printf(1,"C Matrix[0][%d]=%f \n",i,temp);
C.push_back(temp);
}
//D
D=*(*(*doubles_list)[2*columns+1])[0];
vpi_mcd_printf(1,"D Matrix[0][0]=%f \n",D);
if (!error) {
//Initialize array
for (unsigned i=0; i< columns;i++){
zx.push_back(0.0);
zdx.push_back(0.0);
zq.push_back(0.0);
}
}else {
vpi_mcd_printf(1,"$runge_kutta:can not read coefficient.\n");
}
return !error;
}else {
vpi_mcd_printf(1,"$runge_kutta:format of coefficient file should be (row ==2*coulumn+2).\n");
}
}
}
static vpiHandle get_module_handle(vpiHandle call_handle)
{
vpiHandle obj=call_handle;
while (vpi_get(vpiType, obj) != vpiModule) {
obj = vpi_handle(vpiScope, obj);
}
return obj;
}
void runge_kutta::get_delta_time(vpiHandle call_handle)
{
s_vpi_time now;
now.type=vpiScaledRealTime;
vpi_get_time(0, &now);
double current_time=now.real;
now.type = vpiSimTime;
vpi_get_time(0, &now);
int units = vpi_get(vpiTimeUnit,get_module_handle(call_handle));
if (old_time !=current_time){
double scale=pow((double)10,(double)units);//キャストがないと駄目
zh=(current_time-old_time)*scale;//
old_time=current_time;
}
}
void runge_kutta::do_integ()
{
double dk,dy;
vector <double > inputs;
s_vpi_value value;
value.format=vpiRealVal;
// unsigned no_of_int=rows;
unsigned no_of_int=columns;//Dec.20.2004;
//現在のVerilogのインプットアレーを読み出す。
for (unsigned i=0; i< no_of_int;i++){//Read Current Input Array from Verilog
vpiHandle word_index = vpi_handle_by_index(m_read_item,read_mem_left_address+i);
vpi_get_value(word_index, &value);
double temp=value.value.real;
inputs.push_back(temp);
}
if (zh==0.0) return;//zhが未設定ならなにもしない。
//Matrix 乗算
for (unsigned mmode=1; mmode <=4;mmode++){
for (unsigned ix=0;ix< no_of_int;ix++){//行
double sum=0;
for (unsigned iy=0; iy< no_of_int;iy++) {//列
//double coeff=*(*(*doubles_list)[ix])[iy];
double coeff=A[ix][iy];//Dec.20.2004
sum +=zx[iy]*coeff;
}
double temp=inputs[ix];
//zdx[ix]=sum+inputs[ix];//行
zdx[ix]=sum+B[ix]*inputs[ix];//Dec.20.2004
}
//積分処理
for (unsigned m=0;m<no_of_int ;m++ ){
switch (mmode) {
case 1 :
dk=zdx[m]*zh;
dy=0.5*(dk-2.0*zq[m]);
zx[m]=zx[m]+dy;
zq[m]=zq[m]+3.0*dy-0.5*dk;
break;
case 2 :
dk=zh*zdx[m];
dy=cc1*(dk-zq[m]);
zx[m]=zx[m]+dy;
zq[m]=zq[m]+3.0*dy-cc1*dk;
break;
case 3 :
dk=zh*zdx[m];
dy=cc2*(dk-zq[m]);
zx[m]=zx[m]+dy;
zq[m]=zq[m]+3.0*dy-cc2*dk;
break;
case 4 :
dk=zdx[m]*zh;
dy=(dk-2.0*zq[m])/6.0;
zx[m]=zx[m]+dy;
zq[m]=zq[m]+3.0*dy-dk/2.0;
break;
} /* endswitch */
}//end for
}//end for
//結果をVerilogのアレーに書き込む
//Apr.18.2005
vpiHandle word_index = vpi_handle_by_index(m_write_item,write_mem_left_address+0);//Apr.18.2005 no_of_int-1);
//Write the result to Write Array in Verilog
double temp=0;
for (unsigned i=0; i< no_of_int;i++){//Read Current Input Array from Verilog
//double temp=zx[i];
temp +=zx[i]*C[i]+inputs[i]*D;//Dec.20.2004
}
value.value.real=temp;
vpi_put_value(word_index, &value,0, vpiNoDelay);//Apr.18.2005
/*
for (unsigned i=0; i< no_of_int;i++){//Read Current Input Array from Verilog
vpiHandle word_index = vpi_handle_by_index(m_write_item,write_mem_left_address+i);
//double temp=zx[i];
double temp=zx[i]*C[i]+inputs[i]*D;//Dec.20.2004
value.value.real=temp;
vpi_put_value(word_index, &value,0, vpiNoDelay);
}
*/
}
//Verilog IFに接続されているアレー、その他のチェック
//一度だけ読み出しアレーについては、ハンドル情報をSaveしておく
bool runge_kutta::get_property_on_verilog_if()
{
s_vpi_value value;
vpiHandle left_range;
vpiHandle right_range;
if (m_read_item == 0 || m_write_item==0) {
vpi_mcd_printf(1,"$runge_kutta: invalid parameter.\n");
return false;
}
if (vpi_get(vpiType, m_read_item) != vpiMemory) {
vpi_mcd_printf(1,"$runge_kutta:3rd parameter should be real array.\n");
return false;
}
if (vpi_get(vpiType, m_write_item) != vpiMemory) {
vpi_mcd_printf(1,"$runge_kutta:4th parameter should be real array.\n");
return false;
}
value.format = vpiIntVal;
vpi_get_value(mode_item, &value);
if (value.value.integer !=(int)(INITIAL)){
vpi_mcd_printf(1,"$runge_kutta:First call should be mode0.\n");
return false;
}
//Read Memory Input for runge_kutta
left_range = vpi_handle(vpiLeftRange, m_read_item);
value.format = vpiIntVal;
vpi_get_value(left_range, &value);
read_mem_left_address = value.value.integer;
right_range = vpi_handle(vpiRightRange, m_read_item);
value.format = vpiIntVal;
vpi_get_value(right_range, &value);
read_mem_right_address = value.value.integer;
if (read_mem_left_address > read_mem_right_address){
vpi_mcd_printf(1,"$runge_kutta: memory range must be Left<Right\n");
return false;
}
if ((read_mem_right_address-read_mem_left_address+1) < columns){//Dec.10.2004
vpi_mcd_printf(1,"$runge_kutta: invalid memory range.\n");
return false;
}
//Write Memory Output for runge_kutta
left_range = vpi_handle(vpiLeftRange, m_write_item);
value.format = vpiIntVal;
vpi_get_value(left_range, &value);
write_mem_left_address = value.value.integer;
right_range = vpi_handle(vpiRightRange, m_write_item);
value.format = vpiIntVal;
vpi_get_value(right_range, &value);
write_mem_right_address = value.value.integer;
if (write_mem_left_address > write_mem_right_address){
vpi_mcd_printf(1,"$runge_kutta:memory range should be LEFT<RIGHT.\n");
return false;
}
if ((write_mem_right_address-write_mem_left_address+1) < columns){//Dec.20.2004
vpi_mcd_printf(1,"$runge_kutta:invalid memory range.\n");
return false;
}
vpiHandle words = vpi_iterate(vpiMemoryWord, m_read_item);
vpiHandle item = vpi_scan(words);
unsigned wwid = vpi_get(vpiSize, item);
if (wwid !=64){
vpi_mcd_printf(1,"$runge_kutta:memory should be real type of array.\n");
return false;
}
words = vpi_iterate(vpiMemoryWord, m_write_item);
item = vpi_scan(words);
wwid = vpi_get(vpiSize, item);
if (wwid !=64){
vpi_mcd_printf(1,"$runge_kutta:memory should be real type of array.\n");
return false;
}
//Finally,Initial Check all done.
return true;
}
//毎クロック呼び出される。
static int sys_runge_kutta_calltf(char*name)
{
runge_kutta* local_ptr=0;
s_vpi_value value;
vpiHandle sys = vpi_handle(vpiSysTfCall, 0);
vpiHandle argv = vpi_iterate(vpiArgument, sys);
vpiHandle file_handle=vpi_scan(argv);//First Parameter
//First parameter
if (!file_handle) {
vpi_mcd_printf(1,"$runge_kutta:missing parameter.\n");
//Mar.25.2004 vpi_free_object(argv);
return 0;
}
//ファイル名取得チェック
if (vpi_get(vpiType, file_handle)!=vpiConstant &&
vpi_get(vpiType, file_handle)!=vpiParameter){
vpi_mcd_printf(1,"$runge_kutta:file name should be string type.\n");
vpi_free_object(argv);
return 0;
}
if (vpi_get(vpiConstType,file_handle) != vpiStringConst) {
vpi_mcd_printf(1,"$runge_kutta:file name should be string type.\n");
vpi_free_object(argv);
return 0;
}
//ファイル名取得
value.format = vpiStringVal;
vpi_get_value(file_handle, &value);
string file_name=value.value.str;
//2nd parameter
vpiHandle mode_item = vpi_scan(argv);
if (mode_item ==0){
vpi_mcd_printf(1,"$runge_kutta:missing mode parameter.\n");
//Mar.25.2004vpi_free_object(argv);
return 0;
}
map<string ,runge_kutta* > ::const_iterator im=file_map.find(file_name);
if (im ==file_map.end()){//新しいファイル処理
local_ptr=new runge_kutta();
file_map[file_name]=local_ptr;
bool success=local_ptr->file_read(file_name);
if (success) {
bool success=local_ptr->check_and_dump_double_list();
if (success){
//3rd/4th parameter
local_ptr->m_read_item=vpi_scan(argv);//READ ARRAY
if(! local_ptr->m_read_item) return 0;//argv はdelete は、Verilogでデリート済み
local_ptr->m_write_item=vpi_scan(argv);//WRITE ARRAY
if(! local_ptr->m_write_item) return 0;//argv はdelete は、Verilogでデリート済み
local_ptr->mode_item=mode_item;
if(local_ptr->get_property_on_verilog_if()) local_ptr->go_on();
}
}
}else {//既にMapに登録されている。
local_ptr=file_map[file_name];
}
if(local_ptr->go_()){//初期化でエラーなかったらModeパラメータで制御
local_ptr->get_delta_time(sys);
s_vpi_value value;
value.format = vpiIntVal;//モードパラメータを取得
vpi_get_value(mode_item, &value);
unsigned mode=value.value.integer;
if (mode==0) ;
else if(mode==1) local_ptr->do_integ();
}
vpi_free_object(argv);
return 0;
}
//ユーザシステムタスクを登録
extern void sys_ruge_kutta_register()
{
s_vpi_systf_data tf_data;
tf_data.type = vpiSysTask;
tf_data.tfname = "$runge_kutta";
tf_data.calltf = sys_runge_kutta_calltf;
tf_data.compiletf = 0;
tf_data.sizetf = 0;
tf_data.user_data = 0;
vpi_register_systf(&tf_data);
}
|