Blog as you go: sigma delta DAC

I have some piezo speakers from another project. That one used a bi-level amp to drive them. I figured it would  be fun to try a tri-level drive, using an H bridge allows you to have +, – and 0V across the device. And for fun, why not make it a direct sigma delta encoder?

It’s going to run on a microcontroller (an arduino). It’ll need very precise timings, so I’ll not be using the arduino environment.

Here’s a first pass in C++ for linux:


#include <iostream>
#include <cmath>
#include <utility>

using namespace std;

float signal(int t)
{
	//Roughly 0-1
	return (1 + sin(t/20. - 1))/2.1;
}

float quantize(float f, int levels)
{
	return min(levels-1.f, max(0.f, floor(f*levels)))/(levels-1);
}

int main()
{

	float integral = 0;

	for(int i=0; i < 200; i++)
	{
		float output = quantize(integral, 3);
		float difference = signal(i) - output;
		integral += difference;

		cout << signal(i) << " " << output << endl;
	}

}

And this produces output like this (mildly prettified):
graph.png

Which looks about right. The code will pretty much stink on an arduino since it’s all floating point. It’s easy to convert to integer code though:

#include <iostream>
#include <cmath>
#include <utility>

using namespace std;

uint16_t signal(int t)
{
	//Roughly 0-1
	return 65535 * (1 + sin(t/20. - 1))/2.1;
}

uint16_t quantize(int32_t i)
{
	if(i < 21845)
		return 0;
	else if (i < 43690)
		return 32768;
	else
		return 65535;
}

int main()
{

	int32_t integral = 0;

	for(int i=0; i < 200; i++)
	{
		uint16_t output = quantize(integral);
		float difference = signal(i) - output;
		integral += difference;

		cout << signal(i) << " " << integral << " " << output << endl;
	}

}

I’ve used a uint16_t for the value, which will effectively represent both positive and negative levels with 32768 being the zero point. Note that the error integral must be both signed and wider since the errors can grow beyond the signal range:

graph.png

Now to port to the arduino. So, I’ll get my Makefile from here.

I’m going to pick pins 8 and 9 on the UNO, that is PB0,1 on the chip for my outputs. and here’s how you get 3 way opposed outputs as +/-/0. To demo, I connected a pair of LEDs in parallel but facing the other way:


#include <avr/io.h>
#include <util/delay.h>
int main()
{
	int i=0;

    while(1)
    {
		if(i==0)
		{
			DDRB = 3;
			PORTB = 1;
		}
		else if(i==1)
		{

			DDRB = 3;
			PORTB = 2;
		}
		else
		{
			DDRB=0;
			PORTB = 0;
		}

		i++;
		if(i > 2)
			i=0;

        _delay_ms(200);
    }
}

So I started the port and BOOM! 😦

It stopped working. The reason was simple: the simple makefile takes one object file and converts it to HEX. Since I’m using sin(), we actually need to link the .o and .a into a .elf file, then convert THAT to HEX. The snippet is:

%.hex: %.elf
avr-objcopy -j .text -j .data -O ihex $< $@

prog.elf: delta_sigma.o
avr-gcc $(FLAGS) -o prog.elf delta_sigma.o -lm

Obvious, really, in hindsight…

So, OK, now to convert the modulator code to the arduino. Lots of things went wrong. But first, here’s the code:

#include <math.h>
#include <stdint.h>
#include <avr/io.h>

uint16_t signal(int32_t t)
{
	float u = t / 1024.f;

	//Roughly 0-1
	return 65535 * (1 + sin(2*3.151592f*u))/2.1;

}

uint16_t quantize(int32_t i)
{
	if(i < 21845)
		return 0;
	else if (i < 43690)
		return 32768;
	else
		return 65535;
}

int main()
{

	int32_t integral = 0;

	DDRB|=32;

	for(uint32_t i=0; ; i++)
	{
		uint16_t output = quantize(integral);
		int32_t difference = (int32_t)signal(i) - output;
		integral += difference;

		if(output == 0)
		{
			DDRB=255;
			PORTB=1; //Output 1 0
		}
		else if(output == 65535)
		{
			DDRB =255;
			PORTB=2; //Output 0 1
		}
		else
		{
			DDRB=255;
			PORTB=0; //Output 0 0
		}
	}
}

What didn’t go wrong? Nothing! I wasn’t nearly careful enough with my ints (only 16 bits on AVR), ints of specific width, overflow and that sort of thing. Also, initially, I decided to output a 0 level by tri-stating the two outputs, so they both float to the middleish. Turns out that didn’t work well since they float extremely slowly (not surprising really!). Forcing them both down to 0 worked much better.

After all that, I then connected a simple RC filter across it so you an see the results:

That’s actually a pretty nice sine wave there! It ought to be: there’s really not much room for nonlinearity and other distortions to creep in. I’ve zoomed in a few levels so you can see how it looks in detail.

It is however really really slow. I’m using full floating point, and a transcendental operation every iteration of the sigma delta encoder. That is really slowing down the cycle time since the AVR isn’t very fast. That accidentally solves the other problem which I’ve made no attempt to make sure every path takes the same number of cycles. But that sin() is dominating so heavily that it doesn’t matter.

And that’s it for part 1: a working sigma delta encoder. For part 2, I’ll make it fast enough to generate audio tones which aren’t simply the sigma-delta encoder transitions (I hope).

Oh also here’s tehe obligatory github link.

 

Advertisements

How bad is switch bouncing really?

TL;DR: it’s worse.

ETA: I think every electronic engineer goes through this at one time or another. You have a switch that bounces, so you debounce it. I mean every knows that switches bounce. The code doesn’t work 100%, so you fix all the bugs and it still doesn’t work. So then you break out the scope and hells bells! How can it be that bad?! I mean everyone says it’s bad and you knew it was bad but it’s BAAAAAAD.

So, I recently had a switch which needed debouncing. It happens to be a Marquardt (seems to be a good brand) microswitch. In terms of current it’s massively overspecced (16A, 250V), but it turned out to be ideal for the combination of size, button stroke, force and of course price. I guess that sort of switch is used for mains interlocks or something, but I’m using it to switch a GPIO input to a raspberry pi.

Here’s how it’s connected up:

IMG_20170504_184225.jpg

I’d like to say I connected it to both power rails to show how the different contacts bounce. Actually my brain disconnected and that’s how I originally connected up the switch in my project. It’s a terrible idea because you get bouncing from both sets of contacts for twice the fun. The potential divider pulls the middle up to half the voltage so you can see which contacts are bouncing and when.

On the plus side, it’s great for illustrating bouncing…

Here’s what the setup looks like:

IMG_20170504_185448

I use the mini vice (which I’m quite unreasonably proud of) to close the switch really, really slowly. And this is what pushing the switch looks like:

SCR01

(Note: I can’t figure out how to take a screenie without the annoying “take a screenshot” dialog showing…)

It’s quite amazing. Note the “pre-bounce” where a disconnect happens 225ms before the switching starts. In fairness, I couldn’t reliably reproduce that, but it does happen occasionally. Thar main transition looks suspicious though:

SCR02

That’s semi-horrid. The 0V contact opening is pretty clean. It then takes about 9ms for the + contact to start closing. It bounces all over the place, and even after the main bounce has stopped, it still takes 4ms for the voltage to stabilise. I guess connecting only 0V and a pullup would work much better. And here’s the bounces in detail:

SCR03

Horrible, but par for the course. You can see it’s not settling cleanly on the contact in between the bounces. And here’s a different one for comparison (RTFM n00b, the fine manual says how to sake screenies properly):

SCR05

This is bad in a rich and interesting variety of different ways. Like the vast majority of traces, there’s no awful pre-bounce. However, the disconnect is unclean and takes 10ms before it starts to rise properly. After the bouncing nominally finishes, it still takes a whopping 100ms to really settle.

Amazing!

Here’s some samples to fill you with horror as well as a couple of fairly standard ones. Some of them have a longer timebase to show a 200ms settling time.

Some of them have some really weird behaviour as they open. I haven’t figured why and the speed at which I actuate the switch doesn’t seem to affect things much. The same can’t be said for releasing the switch. The slower you release it, the gungier it is. See the enormous 100ms timebase:

yuck yuck yuck 😦 Faster actuation (and the actual bouncing) looks much as you’d expect:

SCR04

though, the time between the break and make is rather shorter.

Either way, the +V contact (the one connected by closing the switch) seems surprisingly much worse than the 0V one. We can of avoid that contact by connecting the switch as a single throw one with a pullup:

IMG_20170504_193440.jpg

And if I’m quick, I can manage to get an open and close on the screen with a 10ms timebase:

SCR05

The push still has some pre-bounce, but it’s nothing like as bad as from the other contact. For the switch being pushed, there is on bounce whatsoever; it’s astonishingly fast:

SCR06

for the release, there’s a bit of bouncing as expected.

Conclusion: don’t connect up switches in a silly way like I did!

That wasn’t in the (LMC555) datasheet :(

So I was making an one-off circuit to drive some things from an RPi and I needed a level shifter. Turns out that a CMOS 555 (on paper) looks like a pretty good bet if you need an ad-hoc solution with mild performance requirements. The TI LMC555 runs all the way down to 3V or so, and can source 100mA from the output. Being CMOS, the output goes more or less to the rails.

So far so good. The way to set it up is to power it off the high side, wire it in Schmitt trigger configuration (pin 6 to pin 2), and set the control voltage at 2/3 of the lower level. And that works just fine.

One problem, it seems that despite being specced to run off 3V, the current sourcing capability drops drastically under about 6V, to the point where at 5V it will only source about 12mA! That’s something of a pity because I needed those 100mA, or more than 12 at any rate, and annoyingly it doesn’t appear to be mentioned anywhere in the datasheet.

😦

On the other hand, I’m glad I bought those jellybeans a while back. I replaced the 555 with a high side PNP switch (a now discontinued BC638 in a small TO92 from one of those Maplin grab bags)  who’s base is driven from a chunky STP55 (a giant TO220) since the latter switches adequate current at only 1V. The 2N7000 is kind of marginal for getting the high base current required to get a low saturated Vce when driven from  3.3V.

So, mission accomplished, but I’m still annoyed about the serious derating. I’ll make a graph when I figure out how to get my scope to act as a data logger.

EDIT: That was way easier than I thought. You can save traces via rather awkward interface to a USB stick in the front panel. So, I set up a 555 to output high into a 10 Ohm load and cranked the supply voltage by hand, measuring the supply and drop across the load. And here’s the result:

The LMC555 hist 100mA at 12V compared to 3 for the LM555

Graph of output pin current into a 10 Ohm load against supply voltage

Well, turns out the LMC555 has a rather high output resistance. The bipolar LM555 on the other hand is a bit of a beast and will give tons of current if you don’t mind the quite high (over a volt) drop at the output.

Jellybeans

I’m low on parts and it’s time to restock. I’m after jellybean prototyping parts—generic ones without any surprising properties. The criteria are they are (a) reasonably cheap (b) available from RS, (c) through hole, (d) available in sensible quantities.

MOSFETS

High power

STMicroelectronics STP55NF06L. Excellent logic level (50A at 3V, over 1A at 2V, something down to 1V), 55A/60V. (about 20p)

Low power

On Semi 2N7000. Logicish level, (50mA at 3V, not much below that), 200mA, 7p.

No outstanding low switching voltage small MOSFETS seem to exist as jellybean parts.

OP-AMPS

OK so everyone says the 741 is too old school and there are opamps better across the board. They’re also 24p which while cheap is actually more expensive than rather better amps. And there are way cheaper ones too.

Super generic

Taiwan Semiconductor TS358CD, 1MHz mostly because it’s 6p. Worse than the 741 in most ways, but 6p! For 2! That’s 3p each! Also, single supply down to 3V. Stretching the definition of jellybean a bit, since it’s so low end, but they’re still going to be fine for basic stuff and if they’re rubbish, well I didn’t lose much. Also did I mention they’re 6p?

Standard voltage

Texas Instruments TIL081, 3MHz, normal voltage range, FET input (33p). A drop-in replacement for a 741, slightly pricier, but better in every single spec. Similar to the LF411 recommended by H&H, not quite as good apparently, but about a tenth of the price.

RRIO, low voltage

Microchip MCP6002-I/P, 1MHz,   1.8V, 100uA single supply, FET? input (1pA) (20p)

Fastish

Microchip MCP6291, 10MHz, 2.4V single supply,  1mA, RRIO, probably FET input (50pA) (35p).

Voltage Regulators

3.3v for things like bluetooth chips and other MCUs. 5V for obvious things and the 317 for everything else.

STMicroelectronics L78L33ACZ, 100mA (14p)

STMicroelectronics L7805CV-DG, 1.5A (9p)

Texas Instruments LM317KTC, 1.5A (15p)

555s

I’ve never tried any of these particular variety of 555s before, but how bad can they be, eh? 🙂 Also yes, I know that an arduino/MCU can do a better job of a timer, but 555s do a bunch of stuff plus there’s always the “what weird things can I do with a 555” game.

Texas Instruments LMC555CN/NOPB, CMOS, 3MHz (74p, kinda expensive for a jellybean, CMOS and fast!)

Texas Instruments NE555P (22p – that’s more like it! Bog standard 555)

 

louder, Louder, LOUDER! (or: more dead bugging)

Sometimes someone makes a chip to do just what you want.

I’ve recently been needing to generate beeps from a BLE113 module (it’s a CC2541) which runs off a CR2032 coin cell at a nominal 3V, but more like 2 to 2.5 in practice. The speaker of choice is a surface mount piezo sounder which are small (9mm square) and unlike the discs don’t require mounting on a sounding board to get sound out. I’ve not idea if those Murata ones are the best, but it’s a respectable brand and those are the first I found that seemed to meet the spec.

They’re not especially loud, only 65dB at 1.5V pk-pk. The microcontroller I’m using has 4 useful channels on timer 1 for this application, and of course the outputs are totem pole outputs. So, driving it with two PWM channels in opposition is driving it with an H bridge which gives the full 2-3V pk-pk swing (depending on the battery voltage).

That makes it little louder, but not an awful lot. The datasheet says that the sounders can be driven at up to 12V pk-pk without damage. The datasheet however merely notes that it is “probable” that increasing the voltage will increase the volume, which is a bit unhelpful, though it has a graph for one (not the one I want) showing an increase with voltage exactly as you’d expect.The question then is how to generate a higher voltage for the buzzer. I had lots of ideas:

    Boost / switched capacitor converter and another H bridge (impractical–too many components)A miniature transformer (none quite small enough or with the right turns ratio)A miniature autotransformer (closer, but still the same problem)Something cunning with an inductor—some sort of ad-hoc boost thing which generates spikes rather than a square wave. Idea not really fully formed.

None of them are really any good. They’re either require impractically large number of components, components that either don’t exist (or I can’t find) or are vague and ill formed and I don’t have the parts to test the idea and anyway I’d probably end up busting up the chip with voltage spikes.

Fortunately it appears that someone thought of this already. It turns out the PAM8904 already does exactly this. It’s a switched capacitor converter with an H bridge, that takes a digital signal in, precisely for the application of driving piezo sounders from low power microcontrollers. Which is nice.

Except I’m not very trusting, and I’ve no idea if it’s worth the effort. I don’t want to order a circuit board and then fiddle around hand soldering QFNs (I’ve seen it done, I’d rather use a stencil) for a one off test. Like so many chips, it’s QFN only now. So the obvious thing to do is to buy one and deadbug it.

I figured I’d try the nice fine hookup wire I’ve got. The colours make it a bit easier to follow which wire is which. Next time, I’d try the same soldering job with enamelled wire. It’s harder to strip and tin, but the insulation doesn’t get in the way. The key to getting the soldering to work in the end was to tape down the wires with masking tape (3M blue tape) as I went along. Even with that it’s two steps forward, one back as you accidentally desolder wires when trying to attach new ones. Here it is!

IMG_20160728_180707IMG_20160728_185724.

(OK, not as good as this, or this, or this—hey that socket is a really nice idea!)

Spot the schoolboy error? I remembered to check continuity between neighbouring pins, but I forgot to pot it or otherwise protect the wires and so some of them fell off when I tried to change the boost voltage selection. And then another 4 wires fell off when I was taking it out. The connection area is tiny and the solder work is frankly not that good, so the joints are amazingly fragile. It’s what I should have done first time, doubly so because the bits of stiff wire for the breadboard really get in the way.

IMG_20160729_134402IMG_20160729_135338

Well, it seems to operate correctly, but I think I’d do it differently next time. A chip socket or veroboard with .1″ header soldered in is a much better choice than flying wires. Potting makes it as robust, but you have to pot it before you know it works.

It’s always a bit hard to tell volume because ears have a logarithmic response and at 4kHz the sound is quite directional. Nonetheless it’s noticeably louder. Yay 🙂

 

 

Good job, TI, I half mean that

I was a little surprised today when trying to debug a board when one of the output voltages was 4.8 V. The main reason for the surprise is that the power supply is specced at 3.3V. And I didn’t make the supply, Texas Instruments did. In fact it’s one of these.

IMG_20160711_182949

Checking the manual reveals that TI do indeed claim that it’s supposed to output 3.3V. Time to find the regulator! There’s no continuity between USB power and the output so it’s not a short. Poking around on likely looking chips quickly reveals that the regulator is OUCH THAT SUCKER IS BOILING HOT! this one:

IMG_20160711_180224.jpg

And has the markings in very small “PHUI”. I got a far as googling “PHUI v” before it autocompleted to “PHUI voltage regulator”, so I guessed I was on the right track :). Not very further down the track the TI TPS730 datasheet crops up showing it’s a TPS73033, which is a 3.3V LDO regulator rated to 200mA. And did I mention it’s baking? It seems to be fried in a rather unfortunate mode (and just for good measure, the NR pin which ought to be at 1.22V is at 0.14). Also, it turns out that the entire circuit diagram was in the manual, so there was no need for that bit of minor sleuthing.

So why good job TI? Well, the other, much more important chips, such as the micro controller I’m programming are rated to 3.9V absolute maximum and it didn’t die with 4.8V across it. I’m pretty pleased about that, because I can imagine going round a cycle frying many chips before finding out the power supply was defective. :shudder: :(.

Well anyway, it’s fried and I can’t use it. I mean technically I’ve successfully programmed the chip and not fried anything as far as I know, but there’s another as yet untested chip on the board rated to only 4.8V absolute maximum. I can’t trust it, so I can’t use it and so as far as I care it’s broken.

And that means I’m free! This guy has a great philosophy which is that if something is broken then there’s no risk of breaking it so you may as well try to fix it. Fortunately, I have some old boards which I’m not currently using which have ST Microelectronics L78L33 SO-8 voltage regulators on them. They’re not LDO so getting 3.3V from 5V is a bit dubious and actually is not quite within spec, but whatever, both my chip and the one in the programmer (a CC2511) work all the way down to 2V, so I reckon that is won’t matter. Also, it’s only rated for 100mA, not 200mA like the original, but both the chips are low power wireless ones so I doubt the current will go too high even when it’s programming. And besides that won’t be for long.

Time to dead bug it! And pot it in hot melt!

IMG_20160711_182349IMG_20160711_182651

And it works! 😎

The LED perhaps to the surprise of no one is dimmer than before and the output voltage is 3.2V which is in fact well within spec for the 7833.

 

Comb filters

I’m currently making some sensitive analog electronics. It’s going to be used in relatively uncontrolled environments and so bits of mains hum sometimes creeps in. The mains is only approximately sinusoidal and so as well as the fundamental at 50Hz, there’s also significant energy in the harmonics at 100Hz, 150Hz and etc.

Digital filters provide a nice way to remove this with a comb filter, so called since it has a comb of equally spaced notches in frequency, so one filter will remove all the harmonics. The basic filter looks like this:

y[t] = x[y] - x[t-k]

Intuitively if you have a frequency component which is exactly k samples long, then the first and second part of the equation will be perfectly in phase since they’re one complete cycle apart and so by subtracting them, that component is removed. The same reasoning applies if there are two, three or more complete cycles in the k steps, so all those components get removed too.

That also applies if there are exactly 0 cycles too, so this filter removes DC.

If you don’t want to remove DC, you can replace the – with a +. The same logic applies, but to integer + half number of cycles instead of full cycles, so that they’re perfectly out of phase.

The derivation of the frequency response is fairly straightforward. Taking Z transforms gives.

H(z) = \frac{Y(z)}{X(z)} = 1 + z^{-k}

Substituting z = e^{j\Omega} and computing the magnitude |H(z)| gives:

\begin{matrix} |H(z)| & = & \sqrt{(1 + e^{j\Omega k})(1 + e^{-j\Omega k})} \\ & = & \sqrt{2 - 2 \cos \Omega k}\hfill \end{matrix}

This gives nulls at a normalised frequency of \Omega k = 0, 2\pi, 4\pi, .... The actual frequency is given by \Omega = 2 \pi f T where T is the sample period, and that gives a fundamental frequency of \frac{1}{kT}. Just to keep it concrete, I’m going to use a sample frequency of 900Hz, and a lag of k=18 giving a filter at 50Hz. Change it to k=15 for a 60Hz filter.

The result is here:

comb-1

Note, just for fun and verification, as well as the theoretical filter response, I’ve plotted the power spectrum of the time domain filter applied to a signal. To make it look nice, the signal needs to be somewhat spectrally flat, so I made a flat spectrum (with a little added noise in the magnitude) and uniform random phase and took the inverse FFT to get the signal.

And back to the filter. While it works, it’s not really very good. The notch is narrow and it affects almost all of the signal, whether near to the notch or not. The problem is that it’s a first order filter. Want to bet that chucking in some higher order terms makes a better filter? How about:

y[t] = x[y] - \frac{1}{2}x[t-k] - \frac{1}{2}x[t-2k]

And the results are:

comb-2

Well that’s promising :). We’ve got sharper notches and a flatter response in the passband (juuust about). So, what’s really going on? If we take a filter like:

H(z) = a z + bz^{-1} + cz^{-2} + dz^{-3}

and crunch through the magnitude of the response, we get:

|H| = \sqrt{\alpha + \beta \cos \Omega + \gamma \cos 2\Omega + \delta \cos 3\Omega}
where:
\begin{matrix} \alpha & = a^2 &+ & b^2 &+& c^2 &+ &d^2 \\ \beta & = ab &+ & bc &+ & cd \\ \gamma & = ac &+ & bd \\ \delta & = ad \end{matrix}

(the general case follows the obvious pattern above). The result is the square root of a Fourier series (in frequency space). In order to get a filter to do just what we want, all we have to do is figure out what we want, square it, compute it’s Fourier series and, er… figure out how to invert the system of polynomial equations above… do that and by the end we ought to have a custom filter.

Trouble is, that’s rather hard :(.

Fortunately this is yet another case where Downhill Simplex saves a huge amount of work.

Aside: I love this algorithm: it has saved me a huge amount of work over the years! It’s a derivative free, unconstrained optimizer. It won’t win any prizes for speed, it’s sensitive to parameter scales and it can quite easily get stuck off a minimum, but in a lot of cases, you can just chuck in a function and let it rip, without messing around with derivatives. It’s fminunc in Octave (or MATLAB), and DownhillSimplex in TooN.

In order to turn it into an optimization problem, we first design the filter in frequency space, then make an error function where the error is the deviation of some sort between the response with a given set of weights and the designed filter. The optimizer then adjusts the weight to minimize the error.

The obvious formulation is simply least squares. However, we can get more flexible very easily. We can rate some areas as more important than others if we like (for example if performance within the notch is more important than in the passband), and we can use the absolute deviation to some power. As the power gets large, the algorithm will tend to minimize the worst case deviation. Note however that Downhill Simplex starts to peter out in its ability to optimize once powers get too large (around 8 or so). The code is pretty straightforward (only tested in Octave):


%The filter is designed by providing a target function on [0, pi).
%Since the filter has only cosines, this is reflected around 0, then
%repeated ad-nauseum
omega=[0:0.01:pi];

%Let the target function start at zero, then step up to 1, but with the
%step being a fragment of some function to make things smoother. Sharper transistions
%lead to larger Gibbs overshoots
%
%Note notches with a slow start (e.g. raised cosine, start > 0) require high order terms
%as features the width of the start region are required, and that's narrow. Ones that ramp
%linearly are very easy as it's very easy to force the functio to be 0 at 0.
width = 0.3;
start = 0.1;
target = ones(size(omega));
ind = find( (omega < width + start) & (omega > start));
target(ind) = .5 - .5 * cos((omega(ind)-start) / width * pi);    %Raised cosine section
%target(ind) = (omega(ind)-start) / width;                         %Straight line
target(find(omega<=start))=0;

% Set the relative importance for different regions. This allows you to trade
% accurace in one region off against another. In this case heavily weight the zero
% part of the the function (i.e. the notch)
importance = ones(size(omega));
importance(find(omega<=start))=10000000;

%Set the number of terms (this is 1 + the order)
n = 40;

% These are the absolute value of the transfer function for
% a filter with the Z transform of:
%       ___
%       \        -i k
% H(z) = >   w  z
%       /__   i
%        i
%
% with a lag of k.

func = @(x, weights)  (abs(exp(j*x'*[0:length(weights)-1]) * weights'))';
phase = @(x, weights)  (arg(exp(j*x'*[0:length(weights)-1]) * weights'))';

%Starting parameters.
%
%Note, I'm interested in designing filters which have the coefficient of
%z positive and all powers have negative coefficients. This is required for
%filters which have f(0) = 0, and importantly, doubles the frequency of a given
%lag, meaning you double the lag size for the same frequency. This gives more
%flexibility with a low sample rate at the penalty of filtering out DC.
%
%So, choose some starting weights that (emperically) do very roughkly the
%right thing, and scale them so the initial function is about the same scale
%as the target function.
weights = [1 -1./ 2.^[1:n-1] ];   %Decay of weights
weights(2:end) = -weights(2:end) / sum(weights(2:end)); %Notch at zero
weights = weights / max(func(omega, weights)') * max(target); % Scale

%Choose an error function which minimizes the sum of some function of the
%absolute deviation, but weighted by the relative importance.
err = @(w) sum((abs((func(omega, w) - target)) .^6).*importance );

%Do the optimization.
wnew = fminunc(err, weights);

We can now get fancy and design custom high order filters. For example we can widen the notch so it has a notch region, but still has steep sides, and require that the filter has good rejection within the notch region and try to reduce the overall worst case deviation. This would help if the inteference is amplitude modulated (i.e. gets better and worse as the person moves), which would cause spreading of the frequencies.

For example here’s a 39th order, 50 Hz comb filter with a widened notch, just because we can. The first graph shows how well the optimized filter matches the design I’ve given it:

comb-3

And here’s what it looks like over the spectrum:

comb-4

Note the 50dB rejection in the entire notch region, the steep transition and the flat passband.

Just for completeness, the weights are (starting at a lag of 0):


  6.9255e-01
 -5.0611e-01
 -3.1259e-01
 -1.6371e-01
 -5.5549e-02
  1.7559e-02
  6.3697e-02
  8.7947e-02
  9.3523e-02
  8.6793e-02
  7.2796e-02
  5.3985e-02
  3.2915e-02
  1.2862e-02
 -4.8328e-03
 -1.9165e-02
 -2.9259e-02
 -3.4763e-02
 -3.6499e-02
 -3.4953e-02
 -3.0909e-02
 -2.5026e-02
 -1.8267e-02
 -1.1082e-02
 -4.1618e-03
  1.9620e-03
  6.5590e-03
  9.8330e-03
  1.1730e-02
  1.2231e-02
  1.1802e-02
  1.1500e-02
  1.0308e-02
  7.8142e-03
  5.8288e-03
  4.6363e-03
  1.8351e-03
 -8.7649e-04
 -3.1003e-03
 -1.8157e-02

They’re still quite large, showing that 39th order isn’t quite enough for the shape we’ve specified, depending on your level of obsessiveness over flatness and available computing power, but the graphs show it’s not bad.

Of course just because we’re filtering in the digital domain and can have filters with vast orders doesn’t mean that’s not overkill. And being essentially a Fourier domain process, all the usual caveats apply, such as sharp edges giving Gibbs overshoots and (with few exceptions—such as sharp notches), small features require very high order filters.

Nonetheless it works, and I it’s useful tool!

As usual: https://github.com/edrosten/filter_design