# Code for sampling and quantizing a signal running too slow

Below is the code for sampling a sine wave and quantizing it with negative feedback. It is running very slowly unless we reduce the sampling size ns.

ClearAll["Global*"]
u[1] = 0; y[1] = 0; ns = 65536; fs = 100; f0 = 13/ns  fs;
x[n_Integer] := x[n] = Sin[2  \[Pi]  f0  (n - 1)/fs];
u[n_Integer] := u[n] = x[n] - y[n - 1] + u[n - 1];
y[n_Integer] := y[n] = If[u[n] > 0, 1, -1];
Do[y[n], {n, ns}]


But the equivalent code in Matlab almost gives the results instantly.

close all;clear;clc;
ns=65536;fs=100;
t=0:1/Fs:(ns-1)/fs;
f0=13/ns*fs;
x=sin(2*f0*pi*t);
y=zeros(size(t));u=zeros(size(t));
for i = 2:length(t)
u(i)=x(i)-y(i-1)+u(i-1);
if u(i)>0
y(i)=1;
else
y(i)=-1;
end
end


Previously I tried to solve the system with RSolve, DifferenceRoot and RecurrenceTable, but did not find a way of defining the quantization(the if condition) in any of the functions. In continuous time domain, we have WhenEvent(in NDSolve/DSolve etc.) which is capable of such defining, but not in discrete time domain.

An illustration of the sampling system is

• 1. Try u[1] = 0.; y[1] = 0.; 2. Have you read blog.wolfram.com/2011/12/07/… ?(Chinese edition can be found here: tieba.baidu.com/p/2186436530 ) 3. Your Mathematica code and MATLAB code aren't equivalent, x, y, u should be defined as List for a fair comparison. 4. Is this your classmate?: mathematica.stackexchange.com/q/301985/1871 Commented May 11 at 11:55
• 1. That simply works! 2. No, but thanks for the sharing.@xzczd 4. No. And it looks a different question. Commented May 11 at 13:02

All numbers in Matlab are machine precision floating point numbers. In Mathematica you can use both machine precision computation and exact computation (and a third thing that is called arbitrary precision). You have to specify on input which one you want to use. So you should make most your numbers floating point like in the following code snippet. Note the small dots after the numbers.

First@AbsoluteTiming[
ClearAll[x, u, y];
u[1] = 0.;
y[1] = 0.;
ns = 65536;
fs = 100.;
f0 = 13./ns fs;

x[n_Integer] := x[n] = Sin[2 \[Pi] f0 (n - 1)/fs];
u[n_Integer] := u[n] = x[n] - y[n - 1] + u[n - 1];
y[n_Integer] := y[n] = If[u[n] > 0, 1, -1];
Do[y[n], {n, ns}]
]


0.402294

Moreover, in Matlab you use arrays to store your data while you use downvalues of symbols in Mathematica. That is similarly inefficient as using cell arrays in Matlab. Certainly somebody told you that good practice in Matlab is to allocated arrays first and then to modify them. You should do the same in Mathematica, too. For generating the vector x you can even use vectorized operations.

First@AbsoluteTiming[
ClearAll[x2, u2, y2];
ns = 65536;
fs = 100.;
f0 = 13./ns fs;

x2 = Sin[(2 \[Pi] f0/fs) Range[0., (ns - 1)]];
u2 = ConstantArray[0., ns];
y2 = ConstantArray[0., ns];

Do[
u2[[n]] = x2[[n]] - y2[[n - 1]] + u2[[n - 1]];
y2[[n]] = If[u2[[n]] > 0., 1., -1.];
, {n, 2, ns}];
]


0.157038

Now the Do-loop is the bottleneck, but that is also a problem in Matlab as in all interpreted languages. If we want to go faster, then we can use Compile to translate the code into C and then compile it:

cf = Compile[{{x, _Real, 1}},
Module[{u, unew, y, ns},
ns = Length[x];
u = Table[0., {ns}];
y = Table[0., {ns}];

Do[
unew = x[[n]] - y[[n - 1]] + u[[n - 1]];
u[[n]] = unew;
y[[n]] = If[unew > 0., 1., -1.];
, {n, 2, ns}];

{u, y}
],
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];


The timings become so short that I have to use RepeatedTiming to get somewhat accurate measurements:

First@RepeatedTiming[
ns = 65536;
fs = 100.;
f0 = 13./ns fs;
x = Sin[(2 \[Pi] f0/fs) Range[0., (ns - 1)]];
{u, y} = cf[x];
]


0.00120244

Can one go faster? In priciple one can write directly in C a use create LibraryFunction via LibraryLink. That allows one, e.g., to avoid haveing to reallocate u and v in every call. I tried it, but for this problem size the overhead of calling the LibraryFunction overweighs; it just does not get faster. But when I increase ns by a factor of 32, I get another speed up of about 3 compared to the method cf`.