matlab_to_c原版完整文件.docx

《matlab_to_c原版完整文件.docx》由会员分享，可在线阅读，更多相关《matlab_to_c原版完整文件.docx（9页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、The MATLAB to C Conversion Process for PAN_TOMPKIN.M Kerry SchutzMay 2015This document covers what changes were made to an entry on MATLAB Central called pan_tompkin.mfor MATLAB Coder compatibility. It is a very popular entry for QRS detection in ECG signals. implementation-ecg-qrs-detectorRequired

2、Products to run pan_tompkin.m MATLAB Signal Processing ToolboxAdditional Required Products to Generate C Code MATLAB CoderYou can quickly find all changes that were made to the original code by searching for the string “TMW” in pt.m. In addition the generated C code was benchmarked against the origi

3、nal MATLAB code for speed. The generated MEX file shows excellent acceleration. All work was performed in R2015a.1. First, I created a copy of the original MATLAB code, pan_tompkin.m, and called it pt.m. Accordingly, I renamed the function entry point to reflect the new name. I also added%#codegen t

4、o allow the editor to do extra C code compatibility checks. This does not impact the execution of the MATLAB code in simulation.function qrs_amp_raw,qrs_i_raw,delay=pt(ecg,fs,gr) %#codegen2. I created a testbench file, tb.m, to exercise the DUT, pt.m. Always thoroughly test your DUT doing anything r

5、elated to code generation. In this case, the author of pan_tompkin.m did not supply a testbench on MATLAB Central nor any test data files toverify against. Since the focus here is on C code generation, I simply substituted random data for the stimulus. Also the authors algorithm has sample rate depe

6、ndence. In this case, Im using 400 Hz for the sample rate but in general you should test against all relevant rates.% tb.m% simple testbench% create stimulus in = rand(1000,1);% call DUTtic; qrs_amp_raw0,qrs_i_raw0,delay0 = pt(in,400,1); toc;3. When I run the tb.m, the DUT appears to work and runs t

7、o near completion but I see an error near the end as well.Conversion to logical from matlab.graphics.Graphics is not possible. Error in pt (line 434)434: if ax(:)435:linkaxes(ax,x);436:zoom on;438: endThe variable ax is a handle graphics axes object and and does not impact the algorithm. To work-aro

8、und the problem, I just removed visualization code from the DUT since this is code best served by the testbench.4. I ran coder.screener on pt.m. This is a course-grain check for MATLAB Coder compatibility. The screener is looking for constructs and functions that are not supported for C code generat

9、ion. It does not take into account data types and sizes which are also important and as such, the build process will have to find those incompatibilities later as we will show in the steps that follow.coder.screener(-text,pt.m,report.txt); edit report.txtEverything looks good at a high-level.5. I cr

10、eated a MATLAB coder project. coder new pt6. I clicked Next and soon discovered the first set of MATLAB Coder incompatibilities.ERROR: Code generation requires variable ax to be fully defined before subscripting it.This first set of errors are all visualization related. Instead of trying to fix this

11、 error, I chose simply to work-around it.WORKAROUND: Removed all visualization code from the DUT. You can quickly find visualization sections of the code by searching for the string “if gr”.Visualization as a general rule should be placed in the testbench and not in the DUT. As a best practice, you

12、should return variables of interest from the DUT function and post-process and/or visualize them in the testbench. Removing visualization code is not however required.Occasionally having visualization within the DUT is convenient for debugging purposes. See the file pt2.m for an example of how to ke

13、ep the visualization within the DUT and still generate C code.7. I specified the testbench name in the project to be tb.m and clicked Autodefine Input Types. Once it determines the types, click Next.8. Next I clicked Check for Issues.9. MATLAB Coder finds more problems.ERROR: Error using Filter Desi

14、gn function, butter. All inputs must be constant.This time it doesnt like the usage of the function butter from the Signal Processing Toolbox. The usage of butter is fine for simulation but not for code generation. The problem is that the frequency specification is not a constant. You cannot generat

15、e C code responsible for butterworth filter design. In other words you cannot use butter as a poor-mans adaptive filter for C code generation. Instead you must design the filter first, offline, and then use the coefficients, online.WORK-AROUND: Make the inputs to butter all constants. In this case,

16、I set fs = 400 Hz but the approach may vary depending on the application. For instance, you may need to pre-define a matrix of filter coefficients for all possible sample-rates your system will use. Or if your filter is truly adaptive and/or the number of samples rates is too large, then going with

17、an adaptive filter approach is probably best,e.g. LMS.fs = 400; % TMW add: butters inputs must be constants10. ERROR: Size Mismatch. Size mismatch (size 0 x 0 = size 1 x 1). The size to the left is the size of the left-hand side of the assignment. The variable qrs_c was initialized as but later chan

18、ged.grs_c is initialized as empty but then changes size (grows) by conditionally tacking on one more sample to a growing vector.qrs_c =; %amplitude of R qrs_i =; %indexSOLUTION: init qrs_c and qrs_i to a 1 by 0 array and make it variable sized, i.e. allowing it to grow.qrs_c = zeros(1,0); %amplitude

19、 of R % TMW add, was qrs_c = ; qrs_i = zeros(1,0); %amplitude of R % TMW add, was qrs_i = ; coder.varsize(qrs_c,qrs_i); % TMW add11. ERROR: Variable y_i is not fully defined on some execution paths.SOLUTION: Add initializations for y_i outside the loop% Thresholding and online desicion rule x_i = 0;

20、 % TMW addedy_i = 0; % TMW added12. It worked! We have now have pt_mex.meww64 and a large pile of ANSI-C code which should both be in perfect agreement with the pt.m.The generated C code was nearly 4000 lines of code including all dependent C files as compared to less than 300 effective lines of MAT

21、LAB code in pan_tompkin.m. But such a comparison is very misleading since MATLAB calls many underlying files which this analysis didnt count. Note that findpeaks.c is 833 lines of code by itself. This is just one line in pan_tompkin.m.f - J Summ 叩Num 归 o f .c fif己19Num 归 o f .h fif es:21U nes o f co

22、de UnesF l F, le deta巾: 3,944: 5,832F;le NameU nes of Codefindpe啦 c833&ml sorl.c826pt.c812&ml_s&top.c221conv.c1特pt_&mx如 l.c126m寸ypes.h95rtGetlnf.c92d,ff c88pt_emxAPl.c85pt_types.he4rtGetNsN.c63filttil t c扔rt_nont;n;te.c43filt” .c39rt nont;n;t&.h37mu n.c29abs.c24pt_&mxuW.h19pt_&mxAPl.h17powe.C16conv.

23、h15abs.h13fli pud. c13,d ff h12eml_setop.h12eml_sort h12fil t臼 h12findp e啦 h12mesn.h12pt.h12filttiIt h11tl;pud.h11P0“ 石h11pt_，m t, s l, ze h11pt _ te, m ; n s t e .h11rt G e tl n f.h10rt G e tN s N .h8pt _ ，m t, s l, ze c7| pt 扫 m;nate.c6Verification and BenchmarkingI wrote a separate M-file to do v

24、erification and benchmarking, verify_mex.m. I compared pan_tompkin against pt.m. Then, assuming those agree, compare pt.m against pt_mex.mexw64. With any luck, all three should provide the same output given the same input.% create stimulus in = rand(1000,1);% call original M DUTtic; qrs_amp_raw0,qrs

25、_i_raw0,delay0 = pan_tompkin(in,400,1); toc;% test for accuracy and speed% to conduct accurate timing benchmarks you need to do multiple runs% and log only the fastest time. A simple for-loop with auto-overwrite does% this.% call modified M DUT for k = 1:10tic; qrs_amp_raw1,qrs_i_raw1,delay1 = pt(in

26、,400,1); Mtime = toc; % repeat multiple times end% call MEX DUT for k = 1:10tic; qrs_amp_raw2,qrs_i_raw2,delay2 = pt_mex(in,400,1); MEXtime = toc % repeat multiple times end% ACCURACY: compare PAN_TOMPKIN to PTisequal_qrs_amp_raw0 = isequal(qrs_amp_raw0,qrs_amp_raw1) isequal_qrs_i_raw0 = isequal(qrs

27、_i_raw0,qrs_i_raw1) isequal_delay0 = isequal(delay0,delay1)% ACCURACY: compare PT to PT_MEXisequal_qrs_amp_raw1 = isequal(qrs_amp_raw2,qrs_amp_raw1) isequal_qrs_i_raw1 = isequal(qrs_i_raw2,qrs_i_raw1) isequal_delay1 = isequal(delay2,delay1)% SPEEDMEX_SPEEDUP_FACTOR = Mtime/MEXtimeThe last section of verify_mex.m benchmarks the speed of pt.m against pt_mex.mexw64. The results show excellent MEX acceleration as one would hope.MATLAB code: 3e-3 seconds per callMEX:3e-4 seconds per call, 10x faster