How can I convert an audio PDM bitstream into analog form?

Question

I have collected PDM bitstreams from MEMS microphones and I'm looking to convert them into analogue form. Does anyone know some practical implementation details on how to do this (beyond a quote from Wikipedia that states this only requires a low pass filter)?

Edit: copied from the comments below

I have some sample PDM bitstreams stored in Excel and I'm looking to convert to amplitude in Matlab.

The row just has an numerical sequence identifier (0, 1, 2, 3, 4, 5, 6 etc...) and then a 0 or 1 for the L and R. For the low pass filtering, how do I pick the design criteria (like the different frequencies)?

The sampling frequency is 4MHz, the max frequency of the sound that can be picked up is 50kHz (mics are being used in an ultrasonic application), PDM bitstreams was captured using digital inputs on a Picoscope oscilloscope. PDM data was not filtered, the raw PDM data from the mics was captured.

Edit: Raw FFT of PDM data:

Answer 1

Since the question contains too few details, this answer is simply a generic filtering code which the OP can use as a starting point. I am posting this since searching below keywords in EE.SE didn't turn up good results for me.

Important functions

1. cheby2 for getting filter coefficients.
2. butter for getting filter coefficients.
3. filter to filter.
4. filtfilt to filter, but adjust for the delay introduced by the filter.
5. lsim can be used for filtering since linear filters are linear systems.

_{Below code doesn't follow good coding guidelines. Use it only as a starting point.}

% expected data format
% column 1 time in seconds
% column 2 0/1 PDM data
data = xlsread('data.xlsx');

t = data(:, 1);
data = data(:, 2);

% find sampling time period of data.
tslist = diff(t);

% check if all the samples are evenly spaced in time
% filtering functions usually expect it to be so.
assert(all(abs(tslist - tslist(1)) < 100*eps));

ts = tslist(1);
fs = 1/ts;

% cut off frequency of low pass filter
fcutHz = 2;
% normalise frequency to lie within 0 to 1
% the filter design function may be expecting it.
fdignorm = fcutHz/(fs/2);

%[numer, denom] = butter(7, fdignorm);
[numer, denom] = cheby2(5, 35, fdignorm);

% filtering using filtfilt
filtered = filtfilt(numer, denom, data);

% filtering using lsim
sys1 = tf(numer, denom, ts);
filtered_lsim = lsim(sys1, data, t);

figure(1);
bode(sys1);

figure(2);
plot(t, [filtered, filtered_lsim], 'linewidth', 1.2);
legend('filtfilt result', 'lsim result');

edit (for more details)

Generating a sample PDM signal

Since a sample of OP's original signal is not available, I have generated a PDM signal in Matlab/Octave which is used in later sections. I used an algorithm from Wikipedia as such for PDM generation. Code below

%% PDM generation

% 4MHz sampling frequency as given in question
fs = 4e6;
ts = 1/fs;

% 50 kHz signal frequency as given in question
f50k = 50e3;
t50k = 1/f50k;

% Let's generate data for 20 cycles of 50kHz
t = [0 : ts : 20*t50k]';

% original signal : 50kHz modulated by a pulse.
os = 0.5 * sin(f50k*2*pi*t) .* (t >= 5*t50k & t<= 15*t50k);

% PDM generation as Given in Wikipedia
% https://en.wikipedia.org/wiki/Pulse-density_modulation#Algorithm
pdm = zeros(length(os));
qe = 0;
for ii = 1 : length(os)
    if(os(ii) >= qe)
        pdm(ii) = 1;
    else
        pdm(ii) = -1;
    end
    qe = pdm(ii) - os(ii) + qe;
end

% change all the -1 to 0 to match format in the question.
pdm(pdm < 0) = 0;

The resulting time domain and frequency domain plots are shown below. Frequency domain plot is taken using abs( fft(.) ).

Analysis of the PDM spectrum before filter design start

We now assume we got the above PDM signal from OP's hardware setup, with only accompanying info being the 4MHz sampling frequency and that the original message is about 50kHz.

We first look at the spectrum of the PDM signal and notice that

1. There is a 50kHz peak corresponding to the original signal and
2. There are some "side bands" due to variations in the amplitude of the original 50kHz signal. The "side bands" are important, but we may have to sacrifice some of them during filtering (trade off).
3. High frequency content especially at 2MHz (Nyquist frequency). I am new to PDM, but this is the strategically placed quantisation noise PDM introduces which can be later removed to recover the original signal.
4. Not shown, but there is a DC content in the PDM signal since the OP is using 0/1 levels instead of -1/1 levels as used in Wikipedia. If we use a low pass filter to recover the original signal, what we recover will have a DC content that was not in the original signal. We can use a band pass filter if we know for sure that the original signal didn't have a DC content. We wont worry about that for now. We will stick to an LPF.
5. In a real life, we won't be able to distinguish the "side bands" from the quantisation noise. We will have to make an educated guess as to where our signal ends and the PDM noise starts. In this case an "empty" region between 1e5Hz and 2e5Hz is probably a good guess.

Constructing the filter specifications

1. Let's assume we want to preserve the main 50kHz and say, four "side bands"; i.e. we want to preserve signals up to 75kHz. The final filter should pass them through without much attenuation; say, less than 0.5dB. The corresponding (pass band) specification is 0dB > |G(f)| > -0.5dB, f < 75kHz.

2. Our filter should also prevent the noise from passing through. We can't get 100% removal. We notice that the noise amplitude is about 350 (51dB) at 2MHz and about 100 (40dB) at 1.1 MHz. We want them to be, say, 12dB lower than the signal main peak which is about 100 (40dB). So we will set our our (stop band) specification as |G(f)| < -23dB, f > 1.1MHz. Similarly we could also look at the attenuation required at other frequencies and see which one forms a stricter specification.

3. Let's also decide that we want the filter gain to go from -0.5dB at 75kHz to -23dB by 300kHz where the noise peaks start to pick up. This will be our transition band / roll off specification. The order of the filter is usually determined by this transition band roll off specification.

4. Depending on the filter type we may need to specify more details like pass / stop band ripple (DSP.SE), phase / group delay at certain frequencies etc. We won't do that here. We have the minimum required stuff now.

5. Since we are realising a digital filter and one of our main noise frequencies is right at the Nyquist frequency, we can get a filter with infinite attenuation there. But we won't worry about that either.

Feed the above specifications into the design equations of the selected filter type or some thing like this. If the equations used are for continuous time filters, we will have to convert it to discrete time. Here, I didn't need any such procedure since the cheby2 function I used directly uses the transition/stop band specification and I just had to select the order by trial and error till pass band response also satisfied the requirements.

I used [numer, denom] = cheby2(4, 23, 300e3/2e6); to get a filter with response shown below (x axis is in rad/s). We can see that it satisfies our requirements above.

Filtered Results

The filtering code at the top of this answer gives the recovered signal. I have shifted (-0.5) and scaled (x2) the filtered result for easy comparison to original. It is shown below. Magnified portions show that the noise is still present. If the amount is too much for the given application, repeat the full procedure above with tighter specifications.