DSPRelated.com
Forums

FFT Speed, FIR Output

Started by groger 7 years ago25 replieslatest reply 7 years ago536 views

Hello,

Hoping there's some user experience here on the CMSIS libraries that could provide some assistance. Forgive my ignorance on some of these questions. As you can tell, I'm not well versed on DSP processing, and I have not written any processing software, but using what's available to me now (CMSIS)

I'm trying to "detect" two short burst sine signal (about 500ms) at 100Hz/150Hz in a noisy environment (random bursts of noise). They are separated by 200ms of time. The signal is input through a low pass anti-aliasing filter, with Fc = 250Hz. Then it goes to an ADC sampling at 1000 samples/sec. On a CMSIS example I found (and there's few to be found), it first uses an FIR filter, and then passes it through an FFT. My first question is on the FIR. 

Using ScopeFIR, I created a "Simple PM" low pass filter with an Fc of -60dB at 200Hz using 120 taps. The question, is the output of the FIR array a "reconstruction" of the signal that's passed into it? Is it possible to take this data and then, for example, scale it to drive a DAC to get the "recovered" signals out? Is the PM filter a good choice? Here is the way I have the FIR configured. 

"OutputBuffer" is a float32 array with size = BUFFER_SIZE.

#define BUFFER_SIZE                     (uint32_t)256

#define BLOCK_SIZE                    (uint32_t)32

#define NUM_BLOCKS                 (uint32_t)BUFFER_SIZE/BLOCK_SIZE

for(uint16_t i = 0; i < NUM_BLOCKS; i++ )
{
arm_fir_f32(&S, &ADC_Samples[0] + (i * BLOCK_SIZE), &outputBuffer[0] + (i * BLOCK_SIZE), BLOCK_SIZE);
}

With respect to the FFT, I realize the FFT will output bins based on the FFT size I've chosen, and this seems to work, if I sample at 1 second intervals. Where the problem occurs is if there, for example, 2 occurrences of the different signals separated by the desired interval of 200ms. If this occurs in the same sampling period of 1 second, it's sometimes not detectable, for some unknown reason. Does this mean using a faster FFT rate? It would be ideal to distinctively detect these occurrences. So this is where I'm spinning my wheels. Using the CMSIS library, how can I sample at a faster rate to detect these 2 signals? 

Or should I forget about using the FFT, then how could I use the FIR output array? Again, I'm trying to understand just what the FIR array contains, there is little documentation that explains this clearly. I guess if your background is in DSP, this might not be a problem either. 8)

Thanks for your guidance.





[ - ]
Reply by Fred MarshallMay 14, 2017

Detection is generally about:

1) Matching, more or less, what's known about the signal to be detected to the input signal.

That's where filtering and/or FFT processing come in.

2) Determining if the "match" is good enough.  This is generally about the amplitude(s) coming out of the matching  process.  

You could use a matched filter but I'm not suggesting that as I don't know enough about your signals of interest.  Also, it appears that frequency variations or range may make that a bit harder to do.  But other methods are going to be nearly as good I should think.

Here's an FFT implementation:

Choose a time frame that's  long enough to contain the temporal extent of the signal of interest.  This will determine the FFT output frequency sample interval.  You could take the magnitude or square the output samples.  With this always-positive output, you can compare it with a threshold to "detect".  The threshold might be made adaptive so that it is some measure above the noise level to reduce false detections.  In fact, this is a good idea.  Too high a threshold and there will be missed detections.  Too low a threshold and there will be false detections.

If you know the bandwidth of the busts then you might add FFT adjacent values to encompass the bandwidth of a burst and then overlap these measures by a factor of 2 so that a burst isn't "split" between frequency windows.  That's assuming that the frequency sample interval isn't already close to that value.  With 500msec bursts, there will be roughly 2Hz bandwidth.  So maybe add 4Hz contiguous samples.

Perhaps the more demanding part will be to detect the presence and the end of the bursts since there is more burst energy than not with 500msec / 200msec spans ON to OFF.  

You ask:

"..is the output of the FIR array a "reconstruction" of the signal that's passed into it? Is it possible to take this data and then, for example, scale it to drive a DAC to get the "recovered" signals out?"

Not exactly.  It's better described as a filtered version of the input.  If the design is reasonable and the input signal is reasonable then it will ideally be the SIGNAL plus LESS NOISE.

If you're going to FFT anyway then you can just ignore the higher frequency components and perhaps do away with the FIR filter - unless it somehow improves the FFT by having a smaller dynamic range to deal with.

[ - ]
Reply by grogerMay 14, 2017

Hi, Thank you for the detailed reply, very helpful. So, in essence, the FIR may not be required?

Then the ADC output's array (outputBuffer) could be input directly to the function:

arm_cmplx_mag_f32(outputBuffer, cfftBuffer, FFT_SIZE);

where "outputBuffer" is presently the sample buffer from the FIR filter, and cfftBuffer is only half the size of the outputBuffer.

You said: "Choose a time frame that's  long enough to contain the temporal extent of the signal of interest.  This will determine the FFT output frequency sample interval."

If my sample rate is 1KHz, how do I choose a "time frame" that is not equal to the time frame of the sample rate? Wouldn't the time frame be 1 second, if the sample rate is 1KHz? Is it as simple as driving the FFT every half-second with 500 samples?

[ - ]
Reply by Fred MarshallMay 14, 2017

The time frame is equal to the *number* of samples "N" times the sample rate.

So, sample rate = 1kHz X number of samples of 1,000 would be 1 second time frame.  But, you could choose something else by increasing or decreasing "N".

[ - ]
Reply by grogerMay 14, 2017

Thanks Fred,

That makes sense. Using CMSIS, I found that the relationship between sampling rate changed the bin where the expected signal was. I think there is other threads in this forum that touched on this. So yes, I will try to sample at different intervals/frame times. It might be better going to a DMA based circular buffer than what I'm presently using.

Roger




[ - ]
Reply by Fred MarshallMay 14, 2017

Well the sample interval in frequency is the reciprocal of N X fs.

So, if you change of of them or both, then the sample interval in frequency will change as well.

[ - ]
Reply by SlartibartfastMay 14, 2017

What are you trying to accomplish with the FIR filter?   The output of the FIR will be the input filtered by the FIRs frequency response.   You've indicated a cutoff frequency of 200Hz at -60dB, so one would hope, but you haven't indicated, that the rolloff is sharp enough that it's not attenuating your desired signals at 150Hz.

The cutoff at 200Hz may not be doing much for you since you indicated that there was already an anti-alias filter at 250Hz.

It's also not clear what you're trying to accomplish with the FFT, so it's hard to provide guidance there, either.   Synchronizing processing to pulse detections can be tricky.   If you have the processing bandwidth available you could do a spectrogram or "waterfall" spectrum and see if that helps.   



[ - ]
Reply by grogerMay 14, 2017

Hi,

I'm more interested in what I can accomplish with the FIR. So let's say I move the anti-aliasing hardware filter to a cutoff of 300Hz. Then I generate a set of digital band-pass filter coefficients with a sharp roll-off that has a band width of 80Hz to cover from 85Hz to 165Hz. After passing through the FIR filter, using the settings and array I have provided, what should I expect to see in the output array? 

This is where I am stuck - I do not know what the output contains relative to the input. Are the samples delayed by the FIR filter (due to the number of taps?) Should the magnitude of the bandpass be the same as the input? Are they re-scaled by the FIR? Can I use them directly into a DAC to generate the signals?

Can the answer please try to be as specific as possible to each question?

Your help is appreciated, thanks!

[ - ]
Reply by grogerMay 14, 2017

And sorry, I did not answer your question: ideally, I would like to have the FIR output to "re-create" the signal (as I alluded to in my latest reply) without the "noise" that's in the spectrum. I'm only interested in seeing the 100Hz and 150Hz (burst) signals.

[ - ]
Reply by SlartibartfastMay 14, 2017

The FIR filter is not doing anything magic.   The output will be the input with the frequencies attenuated by the filter suppressed.   That is all.   If the input has little to no noise or energy outside of the passband, then the output will look very much like the input.

Yes, there will be a delay, equal to half the length of the FIR if the coefficients are symmetric about the center.

The DC gain of the FIR will be equal to the sum of the coefficients.   The gain at the passband frequency will depend on the frequency response of the FIR filter.   If the gain at the passband is unity, then the output scaling will be same as the input.   This is often adjusted using a scalar at the output or by scaling the coefficients.

Since the bandwidth of your described filter passes much of the input bandwidth, the noise suppression will be limited.


[ - ]
Reply by SteveSmithMay 14, 2017

So let's say I move the anti-aliasing hardware filter to a cutoff of 300Hz. Then I generate a set of digital band-pass filter coefficients with a sharp roll-off that has a band width of 80Hz to cover from 85Hz to 165Hz. After passing through the FIR filter, using the settings and array I have provided, what should I expect to see in the output array? 

You should expect that the FIR output has the same amplitude as the input, but will be delayed by one half the number of taps (120/2 x 1 msec/sample = 60 msec).  This assumes that the filter design software hasn't done anything funny.  

The simple PM is fine... I like the windowed-sinc, but at the level you are working the choice of filter won't make any difference. 

The number of taps (which corresponds the width of the transition between the stop and pass bands) can be adjusted to optimize the filter, but it won't be critical. Start with about 120 and see if you can get it working.  If you can't, changing the number of taps is unlikely to make it work. 

Yes, you can run the FIR filter output into a DAC, and it should be just as you expect. However, for a signal at 150 Hz, and sampling at 1,000 Hz per second, that is only about 6-7 samples per cycle.  This may produce too rough of output for what you want. There are several ways to correct this in software, but the simplest way is to increase the sampling rate, say to 10 kHz.  This increases the computation time of the filter, but you probably have plenty to spare.   

Will FIR filtering work?  I assume you have looked at the signals on an oscilloscope and the signal is buried in the noise.  The FIR filter will pass the signals unaltered. Assuming the noise is evenly spread between DC and 300 Hz, retaining the band between DC and 85 Hz will result in about 28% of the signal power getting through, which is about 50% of the noise amplitude. So you should only expect about a factor of 2 improvement in signal-to-noise ratio. Is this enough to help you?  

Regards,

Steve


[ - ]
Reply by grogerMay 14, 2017

Thanks for the detailed reply Steve. I'm going to try some of these things, including the filter type.

There's little to no guidance or good examples (documented or code) using CMSIS libs, so I'm not sure if they expect using this to be intuitive, obvious, subverted, or maybe not for people like me. Anyhow, I am chewing my way through, thanks to help from this forum. Appreciated!

Roger

[ - ]
Reply by CedronMay 14, 2017
Hi Roger,

As I read it, it is not clear from your post whether you are trying to implement a tone detection solution or just using that as an example to learn how the FIR function and DFT work.

If it is the latter, I recommend that you plot your signal and the output of the FIR so you can visually see the results of the operation.

If it is the former, I have two suggested approaches for you.

The first I have implemented, and I know it works.  In your example use a DFT frame size of 32.  For your 100 Hz signal sampled at 1kHz you have 10 samples per cycle.  For your 150Hz signal, you have about 7.  In a 32 bin DFT, these correspond roughly to 3 or 5 cycles per frame.  Thus you only need to calculate a subset of bins, say 3, 4 and 5.  Apply equation (19) from my blog article titled "Exact Frequency Formula for a Pure Real Tone in a DFT":

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

If the result is a real number, meaning the complex component is small compared to the real component, check the real part to make sure it is between -1 and 1 and if it is take the inverse cosine.  The resulting value called $\alpha$ is the frequency in radians per sample.  You can convert this to Hz using your sampling frequency.  If it is close to your expected value then you have a valid tone for that interval.

The second approach is much more computationally efficient, and conceptually understandable.  However, I haven't actually tried it because I just derived the equations within the last few weeks.  The first set of equations (there are more coming) are written up in my latest blog article I posted just today titled "Exact Near Instantaneous Frequency Formulas Best at Peaks (Part 1)"

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

These are formulas that calculate the frequency of a pure tone in the time domain with a "built in FIR".  However, I have not tried them yet in this context, I've only done a bit of testing.  You could be the first to put them to practical use.

There is no need to process your signal in blocks.  Start at the beginning, plus your k choice, and run to the end, minus your k choice.  At each local maximum you find, apply the frequency formula.  If the quotient is between -1 and 1, take the inverse cosine.  The result is also $\alpha$, just like the DFT case.  If the frequency is in you range you know you have a tone in that neighborhood.  With this approach you are getting tone detection resolution down to the cycle size level.

Since this is a learning exercise for you, I've tried to keep my tactical recommendations to a minimum.  I've given you the strategy.  If you have any questions, you can post them in this thread or email me at cedron at exede dot net.

Hope this helps,
Ced
[ - ]
Reply by grogerMay 14, 2017

This is excellent, I will read your blogs and keep you posted!

Roger


[ - ]
Reply by dgshaw6May 14, 2017

Hi Roger,

Do you really need to reconstruct the signal or is it really a detection problem?

My best suggestion for your detection problem is to use a Goertzel detector for each of the tones.  It is a VERY efficient implementation and very accurate discriminator.  Basically a couple of multiplies and a few adds.

Lots of the discussion here has been about various numbers of bins for the FFT.  A Goertzel detector is effectively a single bin detection for an extremely accurate Fourier estimate of the desired frequency.  The accuracy of the selected frequency is constrained only by the quantization of one coefficient - 2*cos(2*pi*desired/sampling).

Goertzel detection is implemented simply as an IIR biquad section.  I like to view it as a sympathetic oscillator.  The poles are on the unit circle and the energy in the internal variables will grow without bound if the desired frequency is present.  The detection process is implemented as an integrate and dump technique.  At regular intervals the detection is run by measuring the magnitude of the integrated content and then resetting these variables and restarting.  The precision of the frequency is defined by the coefficient mentioned above, and the length of the integration process.

I have used these Goertzel detectors for years for DTMF (tone dialing) and call progress tone detection.  Basically 8 DTMF tones in two frequency ranges are used to represent one of 16 possible combinations.  The minimum on time is 60 msec and min detection is required for tones longer than 40 msec.  Call progress is a limited combinations of 350, 440, 480 and 620 Hz.

I hope this helps even if it has nothing to do with the libraries you are looking at.

Regards,

David

[ - ]
Reply by grogerMay 14, 2017

David,

You are correct, I do not need to reconstruct. I was saying that only in the context of trying to get a grasp on the how, and how-to of the CMSIS libs, and in general, the FIR. 

It's is about detecting. Your explanation is very useful as to how the algorithm works. Can you provide some insight from your experience as to how the Goertzel detector works in a  noisy environment, with random adjacent signal frequencies contributing to the noise?

Thanks!


[ - ]
Reply by dgshaw6May 14, 2017

Hi Roger,

The Goertzel detector is the most effective solution that I know of with respect to noise.

The items that you really need to look at are essentially specs.  

In the case of the expected tones, you need to know the exactness of their frequencies and their expected length in time.  You have said about 500 msec which is fine, but have not specified the accuracy of the tone generation process or the accuracy of the A/D converter clock in both PPM and jitter.  All these contribute to the performance, but have similar effects no matter which of the many options you may choose.

You need to know if there are other close frequencies that could cause false triggering, and how long they will last if they occur.

The integration time defines the precision of the detector.  Short obviously means less accuracy, while long means tighter frequency detection, but more ambiguity in the time envelope.

What often happens though is that you have to compromise.  I have had to run multiple detectors in parallel with staggered start times, and then to watch how many trigger in sequence to decide if the tone lasted the desired detection length.

If the precision of the generators are somewhat sloppy, then I have had to run several detectors (3, 5, 7 ...) with frequencies at very close spacing.  This gives you more detection bandwidth without necessarily false triggering on close frequency interference.

I have other ideas, if these are not sufficient to get you started.

Regards,

David

[ - ]
Reply by Rick LyonsMay 14, 2017

Hello groger. I'm assuming you have a digital signal sequence comprising noise-only samples and then suddenly the signal is a 100 Hz (or 150 Hz) sinusoid for roughly 500 msec, and then the signal returns to being noise-only samples. If you want to "detect" when that sinusoidal signal "burst" occurs I suggest you forget about using an FFT algorithm. Instead, I suggest you pass your digital signal through the "envelope detector" shown in Figure 2 described in the blog at:

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

(Although the signals in that Figure 2 appear to be analog signals, they're actually digital samples.) In your case, the output of that Figure 2 envelope detector will be a signal whose amplitude is low in the absence of an input sinusoid, but the detector's output amplitude will suddenly go "high" when an sinusoid is present at the detector's input. When the input sinusoid disappears (goes away) the detector's output values will return to its low-amplitude value.  So, when a burst of sinusoidal signal occurs for 500 msec the detector's output will be a roughly 500-msec pulse (the "envelope" of your sinusoid).  In your case I suggest you start out with the lowpass filter in Figure 2 having a cutoff frequency of, say, 50 Hz.  (You can experiment with different lowpass filter cutoff frequencies once you get the envelope detection process working in your software.)

Once you get the envelope detection process working, you then monitor the amplitude of the detector's output.  When the output amplitude exceeds some amplitude threshold you'll know that a sinusoidal "burst" arrived at the detector's input.  When the detector's output amplitude falls below that amplitude threshold you'll know that the sinusoidal burst is now absent at the detector's input.  If you're able to get a Figure 2 envelope detection process working, you can experiment with some of the other envelope detectors described in the above blog.

Good Luck.


[ - ]
Reply by grogerMay 14, 2017

Hi Rick,

This solution also seems plausible, although I'll likely choke in getting from the math into working code. 8(

Post filter, in addition to going into an ADC input, the signal also pass through a Schmitt trigger shaper to get a square wave, although again, quite a few of those signals are noise. Would it be correct to say the envelope detector (as does the Goertzel algorithm) function "correctly" after the presence of multiple sine signals, and not just random/flicker noise?

BTW - are the envelope detectors in your book? (which I have never gone though completely)

Thanks,

Roger

[ - ]
Reply by CedronMay 14, 2017

Hi Roger,

In light of the presence of the harmonics, you will want to adjust your DFT size so that both your frequencies and their harmonics are whole integer multiples of your frame size.

In 20 samples, the 100Hz signal will have 2 cycles and the 150Hz signal will have 3 cycles.  The harmonics will be aliases and will overlay your bin values so you will want a larger frame.

I would think that 80 samples would be good enough, you'll have to experiment and see.  Anyway, if all the frequencies are integer values in terms of cycles per frame then you won't have any interference from leakage from the harmonics.  Simply look for spikes in the bins of your fundamental.

Ced

[ - ]
Reply by CedronMay 14, 2017

Hi Roger,

Why are you trying to convert the signal to a square wave and then rounding off the corners with a low pass filter?

All you are doing is introducing a bunch of harmonic tones.

This will muddle up both approaches I recommended which are based on the signal being a pure sine wave.

I have a feeling that no matter which detection method you use, the Schmitt trigger shaper will do more harm than good.

To clarify Rick's assumptions, I assumed that you had a relatively clean signal punctuated by noise bursts that could occur anywhere, including within your tone intervals.

Ced

[ - ]
Reply by grogerMay 14, 2017

Hi Ced,

This was an alternate approach using edge detection/capture. This is on a microcontroller with the inputs and abilities to use, so I thought I would try it. That shaping and input are completely separate from the ADC that does any ADC processing.

Yes, assumptions are correct, the noise occurs anywhere in and around the frequencies of interest.

Roger


[ - ]
Reply by CedronMay 14, 2017

Upon reconsideration, I may have been wrong before.  Maybe you do want to keep the squarish waveform and use the harmonics as confirmation.  As I said in my other reply, if you can define your DFT frame so it is an even multiple of your expected wavelengths then the tones and all the harmonic frequencies will align on bins.

In my other post I should have said the higher harmonics will be aliased.  Even with a 20 bin DFT the lower harmonics will still be below the Nyquist frequency.  So without the square wave maybe you can be that granular.

Ced
[ - ]
Reply by Rick LyonsMay 14, 2017

Hi.  The mathematics of the Figure 2 envelope detector that I referred to is very simple.  But you do need to know: (i) what is a lowpass digital filter, (ii) how lowpass filters work, (iii) how to design a lowpass filter, and (iv) how to implement a lowpass digital filter. 

And "No", the envelope detectors in the blog I referred you to are NOT described in my "Understanding DSP" book. (What little I know about envelope detectors I learned after I wrote my "Understanding DSP" book.)

[ - ]
Reply by Rick LyonsMay 14, 2017

Roger, matched filters, Goertzel, FFTs, envelope detectors, etc., Sheece!  People are making all sorts of suggestions to you based on what they think you're trying to do.  But I don't think you've made it perfectly clear what is the exact nature of your input signal and what you are trying to accomplish.  When you use the word "detect", that word means different things to different people.

I think there's hope that someone here can help you solve your problem, but you have to do three things first: 

(i) Post a moderately-detailed block diagram of your system starting from the original source of an analog signal all the way to the A/D converter (include everything you know about your system!).

(ii) Tell us everything there is to know about your original analog signal(s) applied to the A/D converter and post a time-domain plot of that analog signal(s) if possible;

(iii) Tell us exactly what you mean by "detecting" your signal(s).

Don't force the good people here to make assumptions about your system, your signals, or what you're trying to accomplish. Tell us everything there is to know regarding your problem.  If you can do these things then, to quote Susan B. Anthony,  "Failure is impossible."

[ - ]
Reply by grogerMay 14, 2017

Hi Rick:

You are right about this. I agree that any number of these solutions may work if implemented correctly, but the first step is for me to get a bit of background on the DSP fundamentals. As you pointed out in your second last post in "Need to know", getting those basic tenets understood will go a long way to helping me realize what questions to ask and what solution makes sense. 

The thread started out with asking about using a CMSIS library question, but I think it was clear to all the experienced DSP practitioners that responded here, the answer is essentially knowing the fundamentals first. Then most of these responses will be more meaningful.

Thanks everyone that was so kind as to offer a suggestion/answer and provide references to other interesting and useful articles I need to read. 

Roger.