DSPRelated.com
Forums

Generate 2KHz damped sine wave using DAC

Started by KousikBarathwaj 6 months ago19 replieslatest reply 5 months ago274 views

I would like to generate a 2KHz damped sine wave from the DAC of my microcontroller (STM32). How to impose the exponential decay lasting for 1second on the 2KHz signal. I presume multiplication of the sine wave and the exponential decay function would produce the intended result. But i am facing difficulties in the practical implementation as i cannot use the DMA for this purpose.  

Note: I already have generated the sine wave with lookup table containing 64 points. 

Any help or insights onto this would be appreciated. 

[ - ]
Reply by neiroberMay 17, 2024

Do you want to generate a repeated damped sine or just one?  One way to generate a damped sine is to apply an impulse to a 2nd order IIR filter having the desired damping factor zeta and natural frequency fn.  See my post on this at:

https://www.dsprelated.com/showarticle/1341.php

Figure 2 of the post shows an example of a damped sine.  Note this method does not require any memory, other than the coefficients of the IIR filter.  If you want repeated damped sines, you could apply repeated impulses.

regards,

Neil

[ - ]
Reply by KousikBarathwajMay 30, 2024

Hello Neil,

Thanks for you insights and really appreciate your help. I am new to DSP and it took me a while to understand the intricacies and functioning of IIR filter. I have gone through your article and Matlab example for 2nd order underdamped system. I need some more guidance to realize this practically, 

I tried with the following input to get the waveform as follows.

zeta= 0.0005; % damping ratio

fn= 2000; % Hz natural frequency

fs= 10000; % Hz sample frequency

nz= 2; % number of zeros

[K,b,a]= so_demo(zeta,fn,fs,nz) 

damped sine_61008.png

My assumptions:

1. I have chosen fs to be 5 times of the fn to respect Nyquist Criteria.

2. Total number of output samples depends on the total number of samples in the input impulse. Since my total period of decay is 1 second i choose input impulse length (N) to be same as the sampling frequency. 

3. I presume b0,b1,b2 are the feed forward path coefficients and a0,a1,a2 are the feedback path coefficient. My difference equation would take the generalized form as below. This would require 5 multiplications to compute one data point.

de_42693.png

Queries:

1. Should i need 10K input samples(1 followed by all zero) to get a output response of 10K?

2. How do i shift the impulse response to the positive axis without a DC gain, similar to below image.

expected sine_18625.png

Best Regards,

Kousik

[ - ]
Reply by neiroberMay 30, 2024

Hi Kousik,

1.  Yes, that is correct.  But the output will be non-zero for longer than 10K samples for a floating-point system -- hence the name IIR (infinite impulse response).  Note that to obtain a single output damped sinusoid, you could just set the initial condition of the first register to 1, and not have any input signal.

2.  I'm not sure I understand the waveform you show.  It does not look like a rectified version of a damped sinusoid.  Can you provide a formula for the signal?

One note of caution:  given that zeta is so small, the IIR filter is very close to being unstable.  This means that signal truncation in a fixed-point filter could cause oscillation or the wrong damping.  You would need to use enough bits to prevent this from happening.

-- Neil



[ - ]
Reply by KousikBarathwajMay 30, 2024

Neil,

Thank you for your prompt response.  Please find the matlab code for the signal. I shifted the sine to positive axis by adding 1 before multiplying with exponential decay function. 

A = 0.5;   % amplitude 
b = 0.4;   % damping factor 
f = 2;     % frequency 
phi = pi/4; % phase shift
t = 0:0.01:10; % time variable 

y = A * exp(-b*t) .* (sin(2*pi*f*t + phi)+1); 

plot(t,y)
xlabel('Time') 
ylabel('Amplitude') 
[ - ]
Reply by neiroberMay 30, 2024

Kousik,

What you describe is not a damped sinusoid in the usual sense.  If we were to persist in using IIR filters, we could perhaps do the following:


1.  Generate a 2 KHz sinewave and add 1 to it.  Call this unipolar signal u

2.  Generate a decreasing exponential signal.  Call it x

3.  Multiply each sample of u by each sample of x.  


You can generate the sine in a number of ways.  As a point of interest, if the sine frequency f0 = fs/4, then the sine values are simply:

sine = [-1 0 1 0 -1 0 ...].

To generate the decreasing exponential, you can use a first-order IIR filter.  Here is an example:

% first-order IIR filter

a1= -.9993;

b0 = 1;

a= [1 a1];

% impulse response

N= 4096;

x= [1 zeros(1,N-1)];    % impulse

y= filter(b0,a,x);     % impulse response


regards,

Neil

[ - ]
Reply by timbly5000May 17, 2024

Hi. What is your sample rate and how big will the total 1s signal be? If you can fit this into a RAM buffer then simply pre-compute the whole waveform ahead of time, including the exponential decay. A one-shot DMA can them be used to play the sample.

Depending on the STM family, you may be able to compute the sinewave and the amplitude scaling on the hoof, sample by sample. This will depend on your chosen sample rate, CPU core speed and whether it has hardware floating point (probably single precision only but this is enough).

Without h/w floating point, real-time sinf(x) will likely be too slow so you will need to resort to a pre-computed sine table (which you can loop) and apply just the scaling factor in real time.

Finally, I agree with Neil's post on using a damped 2nd-order IIR; this will require hardly any storage but will again require real-time math.

Hardware floating point will make life easier. If you dont have it you'll need to use int32 or int64 numbers and this will bring you into the exciting world of integer-math DSP. You'll need to choose a scaling factor for each of your numbers to provide 'fixed-point' arithmetic with an appropriate balance of range and resolution. Be careful with any IIR solution when using limited-precision math (this includes single-precision floats)

Good luck.

[ - ]
Reply by KousikBarathwajMay 21, 2024

Hello, 

Thanks for offering to help.  I am using STM32L4 series controller. My sinusoidal frequency is of 2KHZ and i have chosen a 64 point sine wave for decent resolution. They decay would last for 1 second and repeat all over.  I tried calculating the points on the runtime but couldn't get the expected output even with FPU enabled. I am calculating these points on DMA half and full complete interrupt callback.  64 points sin wave is pre-calculated and i am only computing the scaling factor in the callback on the runtime. htim1 is configured to tick with 1ms frequency upto 1 second. 

void HAL_DAC_ConvHalfCpltCallbackCh1(DAC_HandleTypeDef* hdac)

{

        uint16_t timer_value = __HAL_TIM_GET_COUNTER(&htim1);

for(int i=0; i<32; i++)

{

Wave_LUT_B[i] = Wave_LUT_A[i] * exp((-0.006) * (timer_value));

}

}

void HAL_DAC_ConvCpltCallbackCh1(DAC_HandleTypeDef* hdac)

{

uint16_t timer_value = __HAL_TIM_GET_COUNTER(&htim1);

for(int i=32; i<64; i++)

{

 Wave_LUT_B[i] = Wave_LUT_A[i] * exp((-0.006) * (timer_value));

}

}

[ - ]
Reply by Robert WolfeMay 17, 2024

Hello,

Depending on how intensive your decay function is, you could possibly have a dual sine wave table.  While DMA is sending to the DAC from table A, a parallel process is modifying table B with your decay - multiplication of existing points by latest decay multiplier - which is installed as the DMA source for next time (in DMA done, or 1/2 done, interrupt).  While B is going out, A is modified, etc. 

Regards,

Robert

[ - ]
Reply by KousikBarathwajMay 21, 2024

Thanks for your comments. That helped me to get started. I am currently using half and full complete DMA callback to compute and fill other part of the buffer.  But still it seems like i am missing as i am not able to get the intended output. 

void HAL_DAC_ConvHalfCpltCallbackCh1(DAC_HandleTypeDef* hdac)

{

   uint16_t timer_value = __HAL_TIM_GET_COUNTER(&htim1);

for(int i=0; i<32; i++)

{

Wave_LUT_B[i] = Wave_LUT_A[i] * exp((-0.006) * (timer_value));

}

}

void HAL_DAC_ConvCpltCallbackCh1(DAC_HandleTypeDef* hdac)

{

uint16_t timer_value = __HAL_TIM_GET_COUNTER(&htim1);

for(int i=32; i<64; i++)

{

   Wave_LUT_B[i] = Wave_LUT_A[i] * exp((-0.006) * (timer_value));

}

}

[ - ]
Reply by Robert WolfeMay 21, 2024

Hello,


I would initially fill the source buffers with simple, deterministic values, and confirm the DMA/signal path first, before trying to deal with the decay part. 

Robert



[ - ]
Reply by Robert WolfeMay 22, 2024

Post deleted by author

[ - ]
Reply by Robert WolfeMay 22, 2024

I would think half-dma you're just changing the DMA source destination for the next time to A or B.  If it was A, change it to B and vice versa. But then if you're sending out A currently, and next time it will be B, B has to be completely updated before the A is done, i.e. before the complete-DMA for A occurs.

Probably a separate task, or better software interrupt. So when A complete-DMA occurs, kick off task/SWI to update A with decay for next time, while B is being transferred.  And vice versa.

Robert

[ - ]
Reply by njswedeMay 17, 2024

This should be trivial.

1) Chose a sampling frequency. 4 times oversampling (16hHz) should be fine.

2) Upon startup, calculate a sample table with something like sin(wt)*exp(-kt).

3) Use the built in timers to synchronize the sample output. 

4) Send the samples one by one the DAC. If you require a continuous signal, keep repeating the samples in the table inidefintely.

This should be less than 100 lines of code. What is the challenge you're facing? Honestly, this sounds a lot like you're trying to have someone do your homework.

[ - ]
Reply by KousikBarathwajMay 21, 2024

Hello, Sorry for sounding too novice and lacking details. I exactly did what you have suggested before but couldn't get the expected output.  My sinusoidal frequency is 2 kHz and the no of points each cycle is 64(To ensure good resolution). My decay should last for a second.  My trigger frequency is 128kHz and i presume it is computationally intensive to calculate sin(wt)*exp(-kt) on the runtime. Precomputing all these points (128k) also seems to not an effective approach.

[ - ]
Reply by Detlef _AMay 18, 2024

I did the work for you. 

It was big fun.

https://www.dsprelated.com/showarticle/795.php


#include <stdio.h>

#include <stdlib.h>

#include <stdint.h>

#include <math.h>

int main()

{

// damped sine

// AMPLITUDE*exp(DAMPING)*cos(2*pi*SIGNAL_F/SAMPLE_F)

#define AMPLITUDE (1234.0)

    float q1 = AMPLITUDE;

#define DAMPING (-0.01)

    float b=exp(DAMPING)*exp(DAMPING);

    float q0,q2=0.0;

#define SAMPLE_F (40000.0)

#define SIGNAL_F (2000.0)

    float a = 2*exp(DAMPING)*cos(2.0*M_PI*SIGNAL_F/SAMPLE_F);

    for(int k=0;k<1000;k++){

     printf("%d %.3f\r\n",k,q1);

     q0=a*q1-b*q2;

     q2=q1;

     q1=q0;

    }

    return 0;

}

[ - ]
Reply by KousikBarathwajMay 21, 2024

Hello, Thanks for offering to help. I tried but only getting distorted waveform at my output. I presume calculating these points on the go at runtime is computationally intensive to be run at higher frequencies,  Since my sinusoidal frequency is of 2KHZ and i have chosen a 64 point sine wave for decent resolution, my DMA trigger frequency would be (2kHz x 64 = 128KHz). 

[ - ]
Reply by timbly5000May 21, 2024

Further observations by me:

You dont say which STM32L4 and how much RAM (or even FLASH) you have.

64 x oversampling is excessive and is challenging your implementation both in terms of computing the samples on the fly and not allowing for pre-computation due to the total sample size being 128k samples (are u using 8 or 16 bit values for the DAC?). An oversampling factor of 8 + an external analogue reconstruction filter might be better?

Another observation is that I wouldn't involve a timer count in the sample value computation; this is potentially bringing a second time source into play and could lead to glitches/discontinuities. I would suggest keep a running count of the # samples since the waveform began in your math routine and derive time from this.

Tim

[ - ]
Reply by artmezMay 17, 2024

It's not fair to expect someone else to supply all the details to your work. You don't even say which specific STM32 you're using and are much too scant on the details of what you have done (like schematic, tables, pictures, maybe an oscope image). So we have to specullate and fill in the details of what you did? And, the devil in in the details!

Engineering requires work, not just knowledge. You need to partition your "design" into manageable steps and identify what you know and what you need to know (e.g. how many bits are needed in the sineware). It can actually be very difficult for an experienced engineer to "guess" what's on the mind of an inexperienced engineer. That can be more difficult then having the experienced engineer do all the work for a complete and valid work statement.

A good way to learn some of that is to read through a bunch of posts here and elsewhere, that document the journey of others in search of knowledge. Enjoy your journey!

[ - ]
Reply by KousikBarathwajMay 21, 2024

Apologies for the lack of details and appreciate your insights. 

[ - ]
Reply by Robert WolfeMay 23, 2024

Please let us know what solution you eventually came up with, I am curious.

Regards,

Robert