Block decimation in matlab 2

[Decimation]

Hi,

I have to decimate one long data sequence that gives me memory problems,i.e. can't do it in matlab. Is this possible using the decimate function on k, m-length sequences such that the outputs cab be combined to give the same result as the decimation of the original sequence. In my case, I have an array of 33110000 doubles that I intend to break up into 100 sub-segments of size 331100 doubles; decimate each by a factor of 340 and hopefully combine the results.

Any thoughts on this matter would be appreciated.

Thanx
Ken

 

[GregLocock]

Not in the general case, accurately. It may be good enough depending on what exactly you are doing -for instance if you were analysing vibration it'd work.

I'd look at paging the data in, somehow, or doing the whole job in a real language. To be honest I'm surprised Matlab can't handle GB of data elegantly. Ah ok, it would need to page the output to disk as well.

Incidentally the matlab decimation function wouldn't like an r factor of 340, 13 is the recommended maximum, with multiple runs as necessary.

If you post in the matlab forum, point back to this thread rather than starting a cross post





Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376

http://eng-tips.com/market.cfm?

 

[2dye4]

Am i nearly correct in that you want to reduce the number of samples in your array of data with the required precondition that the bandwidth of the original sequence must be truncated.

There is no benefit to separating the sequences only problems.

Decimation can only be performed when there is excess bandwidth OR bandwidth that can be sacrificed.

The general rule is to filter your data with a digital filter with sufficient attenuation and cut off frequency such that you can re sample the result without any significant power beyond the shannon limit for the new sample rate.

33 e6 doubles should be 264 megabytes. Matlab should do this.

there is a confusing line in the help pages of Matlab about decimation.

"y = decimate(x,r,n) uses an order n Chebyshev filter. Orders above 13 are not recommended because of numerical instability. In this case, a warning is displayed."


The 13 Greg refers to only applies to the form above when a 13'th order Chebyshev filter is used.

You can safely decimate more than this in one step provided your filter is well designed.

look here for more

http://www.mathworks.com/help/toolbox/signal/ref/decimate.html

 

[Decimation]

Thanx Greg and 2dye4 for your comments.

Thins is the first time that I have used this forum and I was surprised by the quick responses.

With regards to my problem, I will be more precise in my description of what I'm doing.

I suspect that the memory problems are due to the 2gig of ram that I'm using with on XP system and perhaps using matlab 6??. I'm not sure. Anyway, right now, I'm using a narrowband signal of duration 1 second that is sampled at 33110000 Hz. I had some success (I thiink), in decimating the 100 blocks separately by taking the remainder of, say, block #1, of length 331100, which is 280 when using a decimation factor of 340, concatenating this to the next block, etc, until all blocks were decimated. Is there a problem using this approach in a logical sense? As well, can this approach be extended to a wide band signal where the decimation factor is a function of the cut-off frequency required? I'm doing this using the form of the decimation function that is defined by

Y=decimate(x,r,n,'fir')

Where x= input sequence
r = decimation factor
n = # of taps in FIR filter

Hopefully these details will be useful in addressing this problem.

Thanx

Ken

 

[GregLocock]

As 2dye4 says 33 million samples isn't many, I'll try it on my webbook in Octave tonight.

The problem with your approach is that your fir will behave slightly strangely at the beginning and end of each block, so the signal you are analysing is subtly different to what you would have got by decimating the entire block in one go. I suggest you try it with just two short length blocks and compare it with analysing both as one block. The difference may not matter if you are using real data.



Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376

http://eng-tips.com/market.cfm?

 

[Decimation]

Greg,

I have performed the block decimation with some successful results. I have no problem with the fir effects at the end of each block. I will test it further with varying decimation rates and more than one frequency. The trick here seems to be the concatenation of rem(331100,340) with the next block. Hopefully I will have some success here. It would also free up a lot of memory by doing this type of decimation "on-the-fly". I will attach some results, especially with the more-than-one frequency case when I generate them.

Thanx for your help

 

[2dye4]

Ken you can decimate data streamwise additionally to blockwise.
Decimation is merely filtering a sequence to reduce its bandwidth so that it can be resampled at a lower rate without upsetting Nyquist.

The FIR filter shifts the data results by N/2 samples. In other words the first full output of the filter is found when the data is completely filled over the FIR coefficients with the value at the middle of the FIR becoming the first output. I know that sounds kooky but i can't figure a way to explain it much better.
I think that is what Greg was getting at.

You can use
filtfilt() to compensate for this as it filters the data in both directions to recenter the original origin but there is still energy smeared out at each end.
The output of a block of K samples through a FIR of N coefficients is N + 2K-1 in length.

Lots of times i have memory errors in Matlab do to a command making a matrix out of something i intended to be an array.
Check the dimensions of your variables to see if they are what you expect.

 

[Decimation]

Hi 2dye4,

Thanx for the info. Yes, this is simply the group delay or in my case the phase delay of the signal since it is narrowband. However, I read in some of the matlab docs that this delay is taken care of in the decimate function? Not sure what gives here? It is common knowledge that an fir filter does incorporate some output delay.i.e. for example, consider the moving average filter. If you are averaging, say, every N pts of a k-pt sequence, then the first N-1 output pts are corrupt. I will check this out further but you made a very good point here.

By the way, what do you mean by stream wise decimation?

As well, I can only obtain 8 megs of disk swap space with my XP system. Do you know how I can improve on this if you've seen this before. These memory problems could improve with a newer version of matlab, I suspect.

Cheers and thanx

Ken

.

 

[GregLocock]

When I get around to it I'll post an example of the problem with decimating blocks rather than contiguous data. I'd emphasise that if you an afford small errors then it probably isn't worth worrying about. For instance with vibration data I'd be happy to throw half of it away in some circumstances (just as well as that is what hanning does).



Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376

http://eng-tips.com/market.cfm?

 

[2dye4]

I mean by stream wise that the data can be read and processed completely sequentially, like reading a small portion from disk, processing, then writing the results back to disk using only enough memory to handle the individual segments read from the disk.

Essentially decimation is to low pass filter then resample the results. So this low pass filtering can be applied to huge data sets as long as care is taken to handle the intervals as though they aren't really separate.

If you are decimating by a factor of 340 you must filter the data so that the new nyquist frequency is Fs_old/(2*340).
I would guess that you need to create a low pass filter with a cutoff at Fs_old/(2*800) with a high order. Make sure the attenuation of this filter is -60db or more at Fs_old then just take every 340th sample of the output as your new sequence.

Try to limit your data to a size that doesn't crash matlab and then run and check the data types using <whos> to make sure the dimensions are as you expect.

 

[GregLocock]

Here's a quick demo comparing decimating a full frame of a signal compared with splitting it into halves and decimating it.

The green curve is the difference between the two approaches.

I did choose the original signal with malice aforethought.


Cheers

Greg Locock


New here? Try reading these, they might help FAQ731-376

http://eng-tips.com/market.cfm?

 

[Decimation]

Hi Greg and 2dye4,

Thanx so much for all the wonderful tips that you provided. They were a great help.

Greg, I will try your test code and observe the results.

As for the block decimaton, I will actually scrap it. I now decimate the entire sequence without any problem. I simply changed the virtual memory that I need, after some experimentation, to a value that made the program work...don't quite understand it, but it worked.

Anyway, it's been a pleasure dealing with you guys.

Thanx again

Ken