Model variable/jittered clock to sample a signal in Simulink - how does variable sample times work?

I'm trying to understand the impact of a clock jitter in an embedded system and its effect on an FFT calculation.

Thus, I'm creating a jittered clock in Simulink, which triggers a sample & hold subsystem, which shall create discrete samples (like in an ADC) of a continuous 20 Hz sine input signal.

Due to the jittered nature of the clock, I can't use a fixed sample rate (or would need to choose a way higher sample rate to model the jitter), thus I chose the auto solver and continuous sample rates.

Sample rate legend:

I tried to follow Model Effect of Temperature and Jitter on Crystal Oscillation Frequency and created two options for the clock generation, hit scheduler and the variable pulse generator, but those seem not to make any difference.

In the model, you can choose to enable or disable the jitter. f0 is the nominal clock frequency, which is 1 kHz in this example.

There's also a counter to visualize the number of clock cycles. Assuming a simulation time of 2.048 s, exactly 2048 clock cycles are generated, assuming no jitter, which makes sense given the 1 kHz clock. On my PC, when I enable the jitter, only 1822 clock cycles are present.

Below you can see the input signal (20 Hz sine), the clock, and the sampled signal:

In the detail view you can see the desired unsteady clock and the sampled signal.

While this looks good at the first glance, I was wondering, how many data points per step are really generated? Can the be answered given that SL chooses the variable sample rate? Is there a way to generate exactly one data point/sample per clock trigger?

If the clock was stable, I'd use a rate transition (which is in the model, see the screenshot). But rate transitions can't have variable sample rates, e.g. via an in-port it seems.

My next steps would be to calculate the FFT on those samples. But I'm not sure how to make sure that there's one sample per clock.

In the model you can see various approaches I've tried to use to calculate the FFT, but I'm not sure if any is correct.

Unfortunately, the Spectrum Analyzer as well as the buffer blocks need non-continuous sample rates, thus the rate transition seems to be the only choice? But that implies a constant sample rate, which is not true...

Another trial was to calculate the FFT in a MATLAB function but not sure if that approach is correct.

Ultimately I want to estimate the error of the frequency bin estimation in the FFT given an unstable clock. I.e., when there's no jitter in the clock, the 20 Hz sine should appear perfectly. But what if there's a jitter?

Can the Spectrum Analyzer block used without discrete sample rates?

Are there any concepts I could use I'm not aware of? Right now I'm using normal rising edge triggers, would function-call triggers change something? Again, I'd like to make sure that there's exactly one data point/sample available per clock edge - just like it would be in the embedded device's ADC interrupt.

Is the problem clear? You need any more information?

Thanks for any input,

Jan

The model is attached for your reference.

[Paul]

Hi Jan,

What exactly is the model for the clock jitter?

For example, w/o jitter the time vector for the 2.048 sec would be

Ts = 1/1000;
t = (0:2048)*Ts;

Is the time vector with jitter to be modeled by adding independent, random samples to each element of t? For example, something like

jitter = 1e-6*randn(size(t));  % use normal random number as in the Simulink model
jitter(abs(jitter)>3e-6) = 3e-6*sign(abs(jitter)>3e-6); % clip
tjitter = t + jitter;

Or does each element of the jittered time vector depend on the value of the previous element?

Also, what is the block that's outlined in green with |FFT| and what is the block to the left of it?

[Jan Kappen]

Hi @Paul,

my model for the clock jitter is to vary the period length of the clock by an additive normal random number.

The results are the same no matter if the Variable Pulse Generator or the Hit Scheduler are used.

> Or does each element of the jittered time vector depend on the value of the previous element?

No, it's a simple additive model, exactly like you described.

> Also, what is the block that's outlined in green with |FFT| and what is the block to the left of it?

Thats the Magnitude FFT. I could've used the "normal" FFT block and calculate the magnitude of the complex output - should be the same. The block left of it is a Buffer. It's used to transform a number of samples over time to a block of samples which can be used by the FFT.

The problem is, that the Buffer needs a non-continuous sample time, exactly like the Spectrum Analyzer. I think there are similar things going on inside the Spectrum Analyzer block. The idea is taken from here: Transform Time-Domain Data into Frequency Domain

I think the main question remains: How do I create non-uniformly sampled data points in Simulink :)

Thanks!

Answer 1

[Paul]

Because the Hit Scheduler can accept vector inputs, I suppose the simplest approach would be to create the jittered sample times before the simulation runs and then inject them into the Hit Scheduler at T = 0. However, I can see this approach being potentially problematic if the nominal sampling period is small and there is a need to process many frames of data (I don't know what the maximum buffere size is of the Hit Scheduler block).

If you prefer a dynamic approach that executes during a runtime like the approach in the question, I might have a solution, but thought I'd suggest the simpler option first.

[Jan Kappen]

Hi @Paul

I realized that sometimes the Hit Scheduler does not create events, e.g. there are missing "clock" eges etc, thus I've changed to the Variable Pulse Generator.

If possible I'd like to stick to the dynamic approach, because I am not yet sure what further requirements on the clock variation will occur in the future.

However, I've found that putting all the FFT or Power Spectrum stuff into the triggered subsystem works. Inside the subsystem, there are really only data points at the triggered times, i.e., I have a dynamic sample rate. But when the signals are routed outside the block, they may be subsampled by SL, using a heuristic I don't understand.

Simplified model (getting rid of all the FFT stuff for now):

Subsystem:

You can see that the jitter is on, and I'm expecting 1822 samples to be created inside the sample & hold subsystem (= the number of rising clock edges).

I'm exactly counting 1822 data points inside the triggered sub system using a counter. Writing those samples out to the workspace confirms that:

But once the signal "sampled signal" is routed out of the triggered subsystem, there are 5694 samples (counted 5691 samples in Simulink, 3 are missed due to the signal being zero):

So that means, the same signal can have different number of samples using a variable step solver, depending on which level of subsystem hierachy it's watched at or saved to the workspace. I was not aware of that.

PS, the Spectrum Analyzer does work inside triggered subsystems, but not outside of them if the sample rate is not discrete:

I think for me the solution is to do everything dynamic sample rate related inside the triggered subsystem.

Nevertheless, I'm highly interested in your approach!

Thanks!

[Paul]

I ended up with an approach using the Hit Scheduler, but in other respects is, I think, similar to where you ended up. The basic idea is to schedule the next update at a nominal 1000 Hz rate, where the next update is 1/1000 + jitter from the current 1000 Hz step.

The Sine wave has a period of 1/1000 sec and stimulates the hit scheduler at with that sample time. The Schedule subsystem is atomic with a sample time of 1/1000. It looks like this

The output is the delt to the next, jittered sample time. The Random Number variance is (1e-6)^2.

The Hit Scheduler does what it does, and I was pleasantly surprised to find that the Hit Scheduler output is high at T = 0 (so no special action needed to collect a sample at the start of the simulation, but that sample will always be at T = 0; how does your approach work for collecting the first sample?).

Here is the Collect subsystem.

The Sine Wave source could just as well be an input to this system. Anyway, the Tapped Delay block collects the most recent 2047 samples and the current sample (so 2048 samples altogether as in the Question). The Counter Limited block counts up to 2047, note that it starts at 0 so when it equals 2047 we have 2048 samples collected.

The Process subsystem is triggered when every 2048th sample has been collected, assuming we want to process every 2048 samples without any overlap. All that I have in this subsystem is a To Workspace block, but I suppose it's where your signal processing would go. I'm not familiar with the the DSP System Toolbox blocks, perhaps they could make some of this simpler, or maybe this approach won't work at all. I don't know.

Here's the plot of the differences in the time vector in simout2, logged from the collection subsystem, for a 6 second simulation. 6000 samples were collected with the time tag of the last sample a hair larger than 5.999.

[Jan Kappen]

Dear @Paul, some very nice approaches here, I did not think about the tapped delay block, only about the buffer and memory blocks... I'll try to use Variable Integer Delay because it has a variable delay via in-port (was no requirement so far, but likely the number of samples shall be calculated dynamically).

Could you upload the model by any chance?

It's interesting that you don't have to feed back the En signal of the hit scheduler as shown in Model Effect of Temperature and Jitter on Crystal Oscillation Frequency. I really need to try to understand the block more in detail.

I think the key is to put the sampling stuff into the triggered subsystem and not try to process it outside the block where apparantely SL sample rates can change.

Big thanks again!

[Paul]

I was originally feeding back the output of the Hit Scheduler block to the En signal, but that was causing me complications until I realized it's not necessary based on how I understand the problem at hand. Basically, the time of the next sample does not depend on the time of the current sample, so the feedback is not necessary. At one point I was using the feedback and trying to compensate for the difference between the actual and nominal hit times in order to schedule the next, random hit time. For example, if the second sample is recorded at T = 0.0009, then we have to deal with the difference between 0.001 and 0.0009 in order to schedule the third hit that is supposed to nominally be at T = 0.002, not T = 0.0009 + 0.001. But I then realized that level of complication is not necessary and settled on the current approach.

I did not look too closely at the Model Effect of Temperature and Jitter on Crystal Oscillation Frequency, so can't say what the differences are there.

The Buffer block from the DSP System Toolbox may offer better functionality. I don't use that toolbox so can't say one way or the other.

Not sure what you mean by the sample times changing outside the triggered subsystem(s). Sample times in the model are determined from the block sample times, events, and the solver (maybe other stuff?). I think I agree that the processing has to be done inside the triggered subsystems because those subsystems are essentially operating at the effective sample times, non constant though they may be.