Sounds very familiar. Does this mean we're getting somewhere? One learns best by hands-on, eh?
Now you're doing a complex multiply of complex frequency sequences, which is equivalent to a time domain circular convolution, and seem interested that the results are what we tried to tell you they would be.
Did you compute an ABS after the multiply? And, really, using IMPRODUCT in frequency *should* generate different complex results with different alignments of the reference or of the signal in time. But, the magnitudes or ABS at each frequency should be the same from one trial to the next. Is that what you're seeing?
Fred
Didn't find your answer? Ask the community — no account required.
B
Bret Cahill
IMPRODUCT looked like it might have something to do with complex conjugates.
(a + bi) * (c + di) =3D ac - bd + (ad + bc)i
Most wouldn't call a simulator with contrived noise "hands on" but it will save all kinds of time and money for any project that does get to the real hands on stage. Not only does it eliminate a lot of unforced errors but even more, it improves judgement on what to watch out for.
It's not that hard to research and/or measure what might be the actual numbers so simulators should be used before doing anything.
IMABS after the INV FFT.
It's straightforward:
make a reference pulse -- a polynomial is fast and easy -- with zeros above and below
in the next column move it up or down some cells for the clean signal with a time shift
make noise for the entire sample time and add it to the clean signal
FFT both the noisy signal and ref.
IMPRODUCT the noisy signal and ref.
Inverse FFT
IMABS or IMREAL then ABS
AVERAGE for magnitude.
Now shift the referance in time to be along side the signal and repeat the above for the same magnitude. It's pretty interesting that an entire time period can be covered with what seems to be one operation.
To see the noise reduction ABS the original ref and signal, AVERAGE both and divide the filtered by the unfiltered. In that one case it reduced the error from a SNR=3D2 pure sine wave noise from 26% error to
3%.
PSR might take 20 - 50 times as long to do that.
Next time I'll check.
That is true even for IMABS of the transforms.
Bret Cahill
B
Bret Cahill
is
.
Actually PSR got the error down to 8% in just one cycle.
F
Fred Marshall
OK. So, I'll add some comments:
make a reference pulse -- a polynomial is fast and easy -- with zeros above and below
***I used a 10-period sinusoid with 10 samples per period which filled about 60% of the data-filled part of the arrays.
in the next column move it up or down some cells for the clean signal with a time shift
***Because of what is going to happen next, this won't matter........ (why do I keep repeating myself on this?) Well, of course it affects the position of the peak in the circular convolution but that's all.
make noise for the entire sample time and add it to the clean signal
***OK
***Plot both the signal and the signal plus noise.
FFT both the noisy signal and ref.
***Well, really you want to compute the "time inversed complex conjugate" of the signal and then FFT that. Since it's real, just flip the order of the column to get the time reversal.
***Don't forget to zero pad enough to avoid temporal aliasing at step. I doubled their lengths with zero-valued samples.
IMABS both results and plot.
IMPRODUCT the noisy signal and ref.
***Yep
IMABS the Product and plot.
Inverse FFT
***OK
IMABS or IMREAL then ABS
AVERAGE for magnitude.
***I'm not sure that works very well. Why not peak pick here? The peak should occur when there's the maximum SNR and it occurs when there's the maximum overlap in the noiseless case. There are likely more sophisticated ways I'm sure.
***Why not use the abs of the FFT and peak pick there? Need to watch scale factor. Likely signal dependent....
***Since it's magnitude you want, check out Rick Lyon's recent article on this topic. I cited it in a recent post here. Helps estimate magnitude when using and FFTs and when frequency isn't right on one of the FFT frequency samples.
Note: if you use a square sine burst as a test case signal then it's really easy to see things in all the results. The noiseless case is interesting - yielding a triangular envelope. And that gives some idea of what happens as you add noise, helps in scaling, etc.
I'd rather use Matlab or Scilab for something like this but it was instructive to use Excel. Not too bad a way for relatively small sequences.
Fred
B
Bret Cahill
and also IMPRODUCT the ref by itself as a "control" to see the error.
Inv FFT ref * ref.
and ref * ref.
It's important to keep track of how much of what you may be calling "padding" because phase sensitive rectification (time domain) multiplies the noise by zero in this time period and comes out looking pretty respectable.
One more comment about PSR:
To compare PSR with the FFT filter it should be necessary to take the SQRT of the PSR output, or, to avoid dividing by zero, the SQRT(ABS(PSR)).
This may or may not be true for the FFT filter.
Both the PSR and FFT filters can be compared with simple rectification.
There may be a proof that shows the FFT filter will beat PSR every time.
That was only to show it works.
Always check for errors, how it looks, the SNR, etc.
That didn't seem to do anything except change the sign of the final inv FFT.
Are you using the highest frequency in the ref to determine that time?
I tried about 10% of that time and was off by 33% in one test. I attributed that to the fact that the CORREL of the noise with the signal was over 0.93%.
Always check.
Would that work for any arbitrary wave form?
Dividing the inv FFT ref by the inv FFT [cleaned up] signal at each point in time produces a column of mostly correct answers with some discontinuities.
Taking the MEDIAN of that will sometimes give better results than AVERAGE.
MAX(ABS())
I'll check it out.
Square waves seem to be harder to FFT filter. Maybe discontinuities are harder to transform.
Fortunately square waves don't happen by accident.
nces.
Everyone has Excel. It is easy to Email a file and show them what they can expect. That's the whole point of modeling. No one wants to spend a lot of money and then be standing around saying "duh" completely clueless as to why its not working.
You can even post a noisy signal on newsgroups so anyone can copy and paste into Excel and take a look at it.
FYI, Scilab and Octave are both free and replicate a large portion of MATLAB's functionality. I use and prefer MATLAB, but there are free alternatives. That type of environment is much better suited for analyzing this sort of problem than Excel is.
Jason
B
Bret Cahill
equences.
to
Thanks. I'll check them out. What I really need is some way to get noise and histograms to get a better idea of what to expect.
Bret Cahill
F
Fred Marshall
I just explained why you can't let that happen..... You get temporal overlap or "aliasing". So, it *must* be avoided - not an option really unless you want to create anomalous examples.
The original sequence is easily determined to be:
Normalized to: sample interval = 1 second; sample rate = 1Hz time window = 256 seconds frequency interval = 1/256 Hz.
and, of course, one can rescale time and frequency to suit.
The embedded signal is at 1/256 Hz with an amplitude of 2*9.2294/256 =
0.0721 and is advanced in time 62 seconds. That is a cosine with phase of 87.68 degrees .. even though the ends are chopped off. I have not tried to determine any phase error due to the chopping off - there may be some but I rather doubt it. It has a DC component of 0.00873 which is due to the ends being chopped off asymmetrically and 2nd and 3rd harmonics noticeable due to the same sharp edges.
If the signal is correlated with a reference cosine by multiplying in frequency, the result is a sinewave shaped output with a peak at 62 seocnds which corresponds to the delay and phase mentioned above. The amplitude of the peak is 9.2294 which, divided by 256 and multiplied by
2 is 0.0721 as above.
Did I get it close enough?
Fred
B
Bret Cahill
Unless you know exactly when the signal cycle starts and ends noise for a short time just before and after the signal cycle time period is unavoidable. This can be reduced by increasing the "padding" but it cannot be entirely eliminated.
Even if you know the curves it would still be ambiguous as to what was supposed to be the noise and what was supposed to be the signal.
The correlation is so high, 0.93 so good filtering could not be expected.
Bret Cahill
F
Fred Marshall
??????
?????
***No it does not do that. The padding is just to add samples so the sequences are long enough to be circularly convolved. The alternative is to NOT do circular convolution, just linear temporal convolution. Of course then we'd not be talking about FFTs. If the sequences are zero-padded, the latter is exactly the same as the former ... which is what you want.
***There really IS NO "FFT filter". There is an implementation of filtering (which can just as readily be purely in the time domain) that is *implemented* in the frequency domain in order to drastically reduce the number of computations necessary in many cases. Otherwise, why do it?
***It's the classical definition of a matched filter and, depending on the waveform, may make a difference.
***No. The highest frequency has nothing to do with it. Of course, you have to select the proper sample interval/rate.
I'm just mindful of the formula for how many nonzero samples come out of a convolution. For one array with M samples and another with N samples it is: M + N -1.
If one does it in the frequency domain then the sequences must have the same length "L" so you can compute the product in frequency element for element.
One way to look at it is this: To circularly convolve two sequences you want there to be a point in the rotation where they don't overlap or, at a minimum, where they overlap on but one sample. If the signal is length "N" and the reference is length "M", then The signal sequence has to have M zeros and it ends up length N + M the reference sequence has to have N zeros and it ends up length M + N and then it's OK for them to overlap on just one sample each and this means we may, if we wish, subtract 1 to get M + N -1 for both sequences.
I started with a signal nonzero over M=100 samples and noise nonzero over N=256 samples sorta with the assumption that the signal was transient and buried *somewhere* in that noise - the alignment, etc. which may be of interest. So, that means there will be 355 nonzero output samples. Since we plan to do an Excel power-of-2 length FFT, that goes to 512. So I padded them both with zeros out to length 512. The zeros only make room for the output..... If you don't do this and still multiply in frequency then you get temporal overlap in the result - which is generally viewed as at least unfortunate...
***I think it depends on the objective. In some broad sense you're measuring a correlation coefficient as a function of time. So, what does that mean in your system context and to your objective is the question.
- If your objective is to determine temporal placement then that means finding a properly registered temporal peak I do believe in just about every case.
- If your objective is to find frequency then ditto but in the frequency domain.
- If your objective is to find amplitude then you might look in either domain but need to be calibrated somehow - which is waveform dependent.
***Dividing????
***Sorry, I meant a gated sine burst - thus square or gate envelope not a square wave. i.e. a rectangular window shape.
Fred
F
Fred Marshall
I did not realize that we were playing games with anomalous situations.
I could just as well give you a signal called "11" and ask what its constituents are. Then answer 1,2,3,5.
No "filter" is going to separate these things out:
1) the first term in the "noise" is a sinusoid which is not resolved by the time frame used. So, all it can do is spill into the adjacent bins.
2) the second term in the "noise" is a sinusoid which is chopped off at
1/2 period - so it splatters to the max.
3) the "noise" isn't random noise so matched filtering is not the optimum filter. So, why are we talking about it. "Good filtering could not be expected" indeed! We aren't even talking about what kind of filter might be optimum for this anomalous case.
4) "noise" is generally meant to refer to random noise. This "noise" signal is more like what one would call "interference".
5) The "signal" contains what is nearly the lowest resolvable frequency for the time window - plus some higher frequency terms. So all it can do is spill into adjacent bins, etc.
There's nothing to be learned from this - except that there can be ill-posed questions and likely, trolls.
B
Bret Cahill
l
lly
=3D
ase
-
d
e
he
y
That was the worst case -- something Japanese engineers can be expected to consider when rebuilding their nuclear power plants.
There is a real duality with conservative engineers spending so much time considering the "worst possible scenario." It may prevent them from trying more things that should be tried.
It can turn into something contrary to the "optimistic" component of Warren Buffet's secret to success, if not the "open minded" part.
Basically engineers need to be able to do a lot of cost-benefit risk analysis very quickly and intuition even judgement won't always help.
But if you know you are going to be hung in a fortnight anyway what have you got to lose?
It probably wouldn't appear on a histogram.
.
uld
There may be a need to sample over more time, more cycles. Keep comparing the new data with the old and see if it's going to converge before everyone dies of old age.
You used a lot of cycles with a few samples/cycle. I used one cycle with a lot of samples.
Taking more samples/cycle helps up to a certain limit. I just tried
8, 16 and 32. The frequencies/cell were halved then quartered to get the same noise and same signal for each test.
But if the number of samples/cycles is limited by the system and you have the time to sample more cycles, then more cycles is the way to go. Then there will be little chance of significant correlation between the noise and the signal.
If you know that the noise and signal are "well behaved" maybe some kind of regression fit would provide more points for the FFT.
Is that a common tactic?
Bret Cahill
F
Fred Marshall
One needs to set both the sample rate and the temporal window according to what's expected.
- the sample rate well above the Nyquist limit .. how much above depends on the application. 10 samples per cycle is a lot really.
- the temporal window large enough to resolve the smallest frequency and/or the smallest distance between frequency components that you must resolve. One can go beyond this in particular applications and one can apply "tricks" which seems well beyond this discussion.
- and, to do circular convolution or "frequency domain filtering" one has to zero pad the data sequences to avoid temporal overlap in the result. I did that with your signal and interferences.
Some observations:
- If you line up a sequence which is "the signal" alongside a sequence which is that same signal plus noise - however you define the noise - such that the signal buried in the noise is lined up with the pure signal, then that alignment represents the one position of the signal/replica/reference that will generally correlate best with S+N.
- So that one alignment generates but one output point (multiply and sum) in a convolution and it should result in the highest output.
- It's not so easy to do convolution in Excel in the time domain - so your notion of doing the exact same thing in the frequency domain with Excel is a good one.
- The matched filter impulse or unit sample response is the time reversal of the signal. The method described above already does that because the objective is to temporally align all the signal components. They align when the last part of the signal just arrives in the filter and the filter is filled with the signal with the first part of the signal at the "end" of the filter. Obviously, this makes little or no difference with a tone burst at a single frequency unless there are very few cycles - as in the case of your polynomial "signal".
To visualize that you might want to do this: Draw a FIR filter as a transversal filter with delay blocks. Input on the left, output on the right as usual. Show the coefficient values along the length of the filter on the diagram as vertical bars. This "plot" of the coefficients goes from left to right and is the time-reversal of the unit sample response of the filter. So, the filter structure matches the signal structure and the unit sample response is a time reversed version of the structure looking left to right.
The output a matched filter with the signals and interferences you gave does produce a maximum at t=0. (This done by padding signal and S+N with 256 zeros before the FFTs).
However, the output of the matched filter doesn't necessarily approximate the signal waveshape. A tone burst will become a tone triangle .. and so forth.
If there are interfering sources with energy at the same frequencies as the signal, then separation isn't generally possible.
For periodic signal in random noise (a VERY common case) take a look at "autocorrelation". I think that's what you're looking for IF that's a reasonable model of the situation. But for a periodic signal in periodic interference as you presented, how do you tell the difference?
Fred
B
Bret Cahill
I'll get back to this later.
Bret Cahill
Join the Discussion
Have something to add? Share your thoughts — no account required.
Didn't find your answer?
Ask the community — no account required
Report Content
You are reporting this content to the moderators. They will look at it
ASAP.