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

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:

• When you say "collected", do you mean the samples are stored in a file on your PC? If so, you can convert to amplitude using a digital low pass filter using a software such as matlab or GNU Octave. Please describe the details of your setup for specific answers. Include data format, sampling frequency etc. – AJN Sep 16 '20 at 11:58
• It literally does just need a low pass filter. What are you really asking? How to design a suitable LPF? – user_1818839 Sep 16 '20 at 11:59
• @AJN Yes I have some sample PDM bitstreams stored in Excel and I'm looking to convert to amplitude in Matlab. Do you happen to know of any examples I can run through that demonstrates this? – dg1271 Sep 16 '20 at 12:58
• I'm trying to implement a low pass filter in Matlab using a butterworth filter as per the suggested answer (and another resource I found). However, how do I select my cut-off frequency, my order etc? How do I know when I have optimized my filter to get the best analog conversion I can? – dg1271 Sep 16 '20 at 16:48
• The added image doesn't say "raw", it says "filtered". It doesn't say "fft" either. It says "time domain". – AJN Sep 18 '20 at 16:51

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

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.

• This posting could use some improvement, as a lot of the question is about how to determine a suitable cutoff. Your code just invents an unexplained number "2" which isn't necessarily suitable. It would seem from a subsequent posting they aren't getting a useful result. – Chris Stratton Sep 18 '20 at 15:44
• From the initial comment from OP, I believed they wanted an implementation example. Later OP provided two frequencies (4M and 50k). I didn't think this was sufficient info to go into specifics like order selection and cut off selection (Hence my disclaimer at the top of the answer). I didn't follow the latest comments until you commented on my answer. I will go through the later comments and see if it is enough to add more details to this answer. – AJN Sep 18 '20 at 16:29