]]>

Note: In earlier versions of the 802.11 standard, the (BCC) OFDM example was given in Annex G. However, an error existed (e.g. in IEEE Std 802.11-2007) in that the FCS was incorrect (see link or search online for "IEEE 802.11-09/0042r0").

The WiFi (802.11) standards documents are free to download from the IEEE webpages (see, for example, here). (Sadly, you are required to register first, which I would describe as irritating and pointless).

Nonetheless, I became interested in throwing together a simple 802.11g (OFDM) modulator. A really useful walk-through of how to do this is provided in Annex I.1 of this document.

However, there are a few gaps in the worked example and some details are perhaps easy to miss. Therefore, I thought I wouldshare the MATLAB scriptsI used for filling in the gaps. These are rather verbose, which I find makes them more useful as a reference, but less useful as anything resembling a usable implementation.

Unfortunately, the main script alone vastly exceeds edaboard's miserly character limit, so please download the full code here: wifi.zip.

As a quick taster, here is the start of the code. This excerpt gets you as far as the end of the STF:

Code Matlab M - [expand] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 close all; clear all; clc; % This script works through Annex I.1 of the 802.11-2016 standard. % Make sure we can see the helper subroutines addpath('helpers'); % ======================== % % == I.1.1 Introduction == % % ======================== % % "This example uses the 36 Mb/s data rate..." rate_info_lut = get_rate_info(); rate_info = rate_info_lut(6); % "... and a message of 100 octets" length_bytes = 100; % ======================= % % == I.1.2 The Message == % % ======================= % % The 72-character message a = ['Joy, bright spark of divinity,' char(10) 'Daughter of Elysium,' char(10) 'Fire-insired we trea']; % [Table I-1] The 100-byte PSDU in hex % Comprises: 24-Byte MAC + 72-Byte DATA + 4-Byte CRC. table_I1 = [... '0402002e00'; '6008cd37a6'; '0020d6013c'; 'f1006008ad'; '3baf00004a'; '6f792c2062'; '7269676874'; '2073706172'; '6b206f6620'; '646976696e'; '6974792c0a'; '4461756768'; '746572206f'; '6620456c79'; '7369756d2c'; '0a46697265'; '2d696e7369'; '7265642077'; '6520747265'; '61673321b6']; % Gradually manipulate PSDU into vector of bytes table_I1 = reshape(table_I1.',[],1); psdu_hex = reshape(table_I1, 2, []).'; psdu_dec = hex2dec(psdu_hex); % Separate PSDU into MAC header, message and CRC mac_header_dec = psdu_dec(1:24); message_dec = psdu_dec(25:end-4); crc_dec = psdu_dec(end-3:end); % Sanity check message content message_char = char(message_dec).'; assert(isequal(message_char, a)); % ====================================== % % == I.1.3 Generation of the preamble == % % ====================================== % % ------------------------------ % % -- I.1.3.1 Short sequences -- % % ------------------------------ % % [Table I-2] Frequency domain representation of the short sequences. % Note: We'll use better precision here: sqrt(13/6) instead of 1.472. p = 1+1j; STF_26 = sqrt(13/6)*[0, 0, +p, 0, 0, 0, -p, 0, ... 0, 0, +p, 0, 0, 0, -p, 0, ... 0, 0, -p, 0, 0, 0, +p, 0, ... 0, 0, 0, 0, 0, 0, -p, 0, ... 0, 0, -p, 0, 0, 0, +p, 0, ... 0, 0, +p, 0, 0, 0, +p, 0, ... 0, 0, +p, 0, 0].'; % Zero-pad up to OFDM symbol size table_I2 = [0*(-32:-27).'; STF_26; 0*(27:31).']; % FFT shift to put DC where it belongs and do OFDM IFFT stf_64 = ifft(fftshift(table_I2)); % [Table I-3] One period of IFFT of the short sequences table_I3 = [ 0 0.046 0.046 1 -0.132 0.002 2 -0.013 -0.079 3 0.143 -0.013 4 0.092 0.000 5 0.143 -0.013 6 -0.013 -0.079 7 -0.132 0.002 8 0.046 0.046 9 0.002 -0.132 10 -0.079 -0.013 11 -0.013 0.143 12 0.000 0.092 13 -0.013 0.143 14 -0.079 -0.013 15 0.002 -0.132 16 0.046 0.046 17 -0.132 0.002 18 -0.013 -0.079 19 0.143 -0.013 20 0.092 0.000 21 0.143 -0.013 22 -0.013 -0.079 23 -0.132 0.002 24 0.046 0.046 25 0.002 -0.132 26 -0.079 -0.013 27 -0.013 0.143 28 0.000 0.092 29 -0.013 0.143 30 -0.079 -0.013 31 0.002 -0.132 32 0.046 0.046 33 -0.132 0.002 34 -0.013 -0.079 35 0.143 -0.013 36 0.092 0.000 37 0.143 -0.013 38 -0.013 -0.079 39 -0.132 0.002 40 0.046 0.046 41 0.002 -0.132 42 -0.079 -0.013 43 -0.013 0.143 44 0.000 0.092 45 -0.013 0.143 46 -0.079 -0.013 47 0.002 -0.132 48 0.046 0.046 49 -0.132 0.002 50 -0.013 -0.079 51 0.143 -0.013 52 0.092 0.000 53 0.143 -0.013 54 -0.013 -0.079 55 -0.132 0.002 56 0.046 0.046 57 0.002 -0.132 58 -0.079 -0.013 59 -0.013 0.143 60 0.000 0.092 61 -0.013 0.143 62 -0.079 -0.013 63 0.002 -0.132]; % Discard columns of indices table_I3(:,[1,4,7,10]) = []; % Convert to complex table_I3 = complex(table_I3(:,1:2:end), table_I3(:,2:2:end)); % Stack in order table_I3 = reshape(table_I3.', [], 1); % Sanity check (up to 3 decimal places) table_I3_err = stf_64 - table_I3; table_I3_err = [real(table_I3_err); imag(table_I3_err)]; assert(all(abs(table_I3_err) < 0.001)); % Extend periodically up to 161 samples. N.B. We don't need to do anything % else for STF due to its 16-sample periodicity. stf_161 = repmat(stf_64, ceil(161/64), 1); stf_161 = stf_161(1:161); % Apply window function stf_161(1) = 0.5 * stf_161(1); stf_161(161) = 0.5 * stf_161(161); % [Table I-4] Time domain representation of the short sequence table_I4 = [ 0 0.023 0.023 1 -0.132 0.002 2 -0.013 -0.079 3 0.143 -0.013 4 0.092 0.000 5 0.143 -0.013 6 -0.013 -0.079 7 -0.132 0.002 8 0.046 0.046 9 0.002 -0.132 10 -0.079 -0.013 11 -0.013 0.143 12 0.000 0.092 13 -0.013 0.143 14 -0.079 -0.013 15 0.002 -0.132 16 0.046 0.046 17 -0.132 0.002 18 -0.013 -0.079 19 0.143 -0.013 20 0.092 0.000 21 0.143 -0.013 22 -0.013 -0.079 23 -0.132 0.002 24 0.046 0.046 25 0.002 -0.132 26 -0.079 -0.013 27 -0.013 0.143 28 0.000 0.092 29 -0.013 0.143 30 -0.079 -0.013 31 0.002 -0.132 32 0.046 0.046 33 -0.132 0.002 34 -0.013 -0.079 35 0.143 -0.013 36 0.092 0.000 37 0.143 -0.013 38 -0.013 -0.079 39 -0.132 0.002 40 0.046 0.046 41 0.002 -0.132 42 -0.079 -0.013 43 -0.013 0.143 44 0.000 0.092 45 -0.013 0.143 46 -0.079 -0.013 47 0.002 -0.132 48 0.046 0.046 49 -0.132 0.002 50 -0.013 -0.079 51 0.143 -0.013 52 0.092 0.000 53 0.143 -0.013 54 -0.013 -0.079 55 -0.132 0.002 56 0.046 0.046 57 0.002 -0.132 58 -0.079 -0.013 59 -0.013 0.143 60 0.000 0.092 61 -0.013 0.143 62 -0.079 -0.013 63 0.002 -0.132 64 0.046 0.046 65 -0.132 0.002 66 -0.013 -0.079 67 0.143 -0.013 68 0.092 0.000 69 0.143 -0.013 70 -0.013 -0.079 71 -0.132 0.002 72 0.046 0.046 73 0.002 -0.132 74 -0.079 -0.013 75 -0.013 0.143 76 0.000 0.092 77 -0.013 0.143 78 -0.079 -0.013 79 0.002 -0.132 80 0.046 0.046 81 -0.132 0.002 82 -0.013 -0.079 83 0.143 -0.013 84 0.092 0.000 85 0.143 -0.013 86 -0.013 -0.079 87 -0.132 0.002 88 0.046 0.046 89 0.002 -0.132 90 -0.079 -0.013 91 -0.013 0.143 92 0.000 0.092 93 -0.013 0.143 94 -0.079 -0.013 95 0.002 -0.132 96 0.046 0.046 97 -0.132 0.002 98 -0.013 -0.079 99 0.143 -0.013 100 0.092 0.000 101 0.143 -0.013 102 -0.013 -0.079 103 -0.132 0.002 104 0.046 0.046 105 0.002 -0.132 106 -0.079 -0.013 107 -0.013 0.143 108 0.000 0.092 109 -0.013 0.143 110 -0.079 -0.013 111 0.002 -0.132 112 0.046 0.046 113 -0.132 0.002 114 -0.013 -0.079 115 0.143 -0.013 116 0.092 0.000 117 0.143 -0.013 118 -0.013 -0.079 119 -0.132 0.002 120 0.046 0.046 121 0.002 -0.132 122 -0.079 -0.013 123 -0.013 0.143 124 0.000 0.092 125 -0.013 0.143 126 -0.079 -0.013 127 0.002 -0.132 128 0.046 0.046 129 -0.132 0.002 130 -0.013 -0.079 131 0.143 -0.013 132 0.092 0.000 133 0.143 -0.013 134 -0.013 -0.079 135 -0.132 0.002 136 0.046 0.046 137 0.002 -0.132 138 -0.079 -0.013 139 -0.013 0.143 140 0.000 0.092 141 -0.013 0.143 142 -0.079 -0.013 143 0.002 -0.132 144 0.046 0.046 145 -0.132 0.002 146 -0.013 -0.079 147 0.143 -0.013 148 0.092 0.000 149 0.143 -0.013 150 -0.013 -0.079 151 -0.132 0.002 152 0.046 0.046 153 0.002 -0.132 154 -0.079 -0.013 155 -0.013 0.143 156 0.000 0.092 157 -0.013 0.143 158 -0.079 -0.013 159 0.002 -0.132 160 0.023 0.023 000 000000 00000 000 000000 000000 000 00000 000000]; % Discard columns of indices table_I4(:,[1,4,7,10]) = []; % Convert to complex table_I4 = complex(table_I4(:,1:2:end), table_I4(:,2:2:end)); % Stack in order table_I4 = reshape(table_I4.', [], 1); % Discard zero-padding at the end table_I4 = table_I4(1:161); % Sanity check (up to 3 decimal places) table_I4_err = stf_161 - table_I4; table_I4_err = [real(table_I4_err); imag(table_I4_err)]; assert(all(abs(table_I4_err) < 0.001));

Again, please download the full code (wifi.zip) and feel free to post any questions here.

]]>

In Part I of this blog, we saw that we can construct two partially correlated signals, x_{1}(t) and x_{2}(t), as:

where P_{1}and P_{2}are the desired signal powers andρis their desired cross-correlation. s_{1}(t) and s_{2}(t) are uncorrelated signals with zero mean and unit power.

We will now generalise this two-signal concept to N signals.

Generalising to N Correlated Signals

By inspecting our expression for x_{2}(t), we can see that (loosely speaking) it has been designed in two parts: the s_{1}(t) part which correlates with x_{1}(t) and the s_{2}(t) part which is uncorrelated with x_{1}(t).

So, if we wanted to design a third signal, x_{3}(t), then we would need a third uncorrelated signal, s_{3}(t), in order to freely choose correlations with x_{1}(t) and x_{2}(t). More generally, we expect then^{th}signal, x_{n}(t), to be a function ofnuncorrelated signals, s_{1}(t), s_{2}(t), ... , s_{n}(t). We can write this idea in matrix form:

or, defining a variable for each bracketed term:

whereAis a lower triangular matrix.

So, for example, in the 2-signal case (i.e. N = 2), we showed in Part I of this blog that:

But how do we find the values of the lower triangular matrix,A, in general? The answer is theCholesky decomposition of the desired theoretical covariance matrix.

To clarify, let's say we want our signals to have covariance matrixR_{xx}:

where (.)^{H}= ((.)^{*})^{T}denotes the conjugate transpose. Then, since all covariance matrices are positive semidefinite, there must exist a lower triangular matrix,L, such that:

whereLis a lower triangular matrix. We can rewrite this as:

whereIis the (N x N) identity matrix. So, by comparing this to Equation (18), we can clearly defineA=Lin Equation (17), so that:

There are several methods for obtainingLin practice. I will first give a full Matlab example, before describing how to compute the Cholesky decomposition of a complex matrix using any programming language you want.

Matlab Example

Code Matlab M - [expand] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 close all; clear all; clc; % Number of samples to generate Nsamps = 100000; % Number of signals N = 4; % Set the required second order statistics P1 = 2.3; P2 = 0.75; P3 = 3.4; P4 = 1.23; rho12 = 0.2 - 1j*0.3; rho13 = -0.6 + 1j*0.1; rho14 = -1j*0.4; rho23 = 0.1 + 1j*0.1; rho24 = 0.5; rho34 = -0.3 - 1j*0.1; % 4x4 (theoretical) covariance matrix Rxx = [P1, sqrt(P1*P2)*rho12, sqrt(P1*P3)*rho13, sqrt(P1*P4)*rho14; sqrt(P2*P1)*conj(rho12), P2, sqrt(P2*P3)*rho23, sqrt(P2*P4)*rho24; sqrt(P3*P1)*conj(rho13), sqrt(P3*P2)*conj(rho23), P3, sqrt(P3*P4)*rho34; sqrt(P4*P1)*conj(rho14), sqrt(P4*P2)*conj(rho24), sqrt(P4*P3)*conj(rho34), P4] % Generate zero-mean, unit-power, uncorrelated Gaussian random signals s = randn(N, Nsamps); % Create partially correlated signals L = chol(Rxx,'lower'); X = L*randn(N,Nsamps); % 4x4 (practical) sample covariance matrix Rxx_samp = X*X' / Nsamps

The above code will produce an output similar to the following:

As in the 2-signal case, note that we never expect the sample covariance matrix to exactly match the theoretical covariance matrix. Sample statistics are themselves random variables, having their own probability distributions.Code:Rxx = 2.3000 + 0.0000i 0.2627 - 0.3940i -1.6779 + 0.2796i 0.0000 - 0.6728i 0.2627 + 0.3940i 0.7500 + 0.0000i 0.1597 + 0.1597i 0.4802 + 0.0000i -1.6779 - 0.2796i 0.1597 - 0.1597i 3.4000 + 0.0000i -0.6135 - 0.2045i 0.0000 + 0.6728i 0.4802 + 0.0000i -0.6135 + 0.2045i 1.2300 + 0.0000i Rxx_samp = 2.2921 + 0.0000i 0.2602 - 0.3927i -1.6626 + 0.2777i -0.0055 - 0.6739i 0.2602 + 0.3927i 0.7513 + 0.0000i 0.1649 + 0.1632i 0.4805 - 0.0006i -1.6626 - 0.2777i 0.1649 - 0.1632i 3.3858 + 0.0000i -0.6106 - 0.2066i -0.0055 + 0.6739i 0.4805 + 0.0006i -0.6106 + 0.2066i 1.2347 + 0.0000i

Other Programming Languages

Free source code is provided in 37 languages for computing the Cholesky decomposition of real matrices at this Rosetta Code page: [link].Please notethat many of these implementations need a very minor modification to work with complex-valued matrices.

Specifically, one term needs conjugating. Using the C code example as pseudocode, we need to conceptually modify the following:

Code C++ - [expand] 1 2 3 4 // Replace this operation... s += L[i * n + k] * L[j * n + k]; // ... with this operation s += L[i * n + k] * conj(L[j * n + k]);

]]>

: [link].: This blog is split into two parts because edaboard limits the number of pictures and equations per blog. The second part can be found hereNote

Introduction

It is often useful to be able to generate (real or complex-valued) signals with specificcross-correlations. In the simplest case, this means we want to produce two signals with a specified cross-correlation (and specified powers). Or, in the general case, it means we want to specify an (N x N)covariance matrixand generate N signals which will produce this covariance matrix.

As a particularly useful example of a type of "signal", we will consider Gaussian (a.k.a. Normally distributed) random variables. However, the ideas presented here may be applied for different types of signals.

In this blog, I will explain in detail how to generate partially correlated signals using common mathematical operations. I will show that discrete correlated signals areextremely easy to generate in software, using theCholesky Decomposition(for which free source code is provided in 37 programming languages at this Rosetta Code page). I will show that in Matlab (or GNU Octave) we basically only need thesetwo lines of code:

where Rxx is the (N x N) desired covariance matrix, N is the number of signals and Nsamps is the number of discrete samples to generate.

I will provide full Matlab source codes for the ideas in this blog. I will also describe the ideas in detail, so you can implement them in any way you want.

If you're not interested in the introductory discussion and derivations for the simple 2-signal case, you can skip straight to Part II of this blog, which contains the main results.

The Simple 2-Signal Case

It is easiest to start by discussing the simple 2-signal case, before generalising to N signals. Let us denote the signals we are designing as x_{1}(t) and x_{2}(t) and assume them to be complex-valued (since it is trivial to then convert the result for real-valued signals if required).To keep the discussion simple, we will assume we are interested in zero-mean signals:

where curly E{...} denotes the expectation operator. So, our task is to design x_{1}(t) and x_{2}(t) to have the following second order statistics:

where * denotes the complex conjugate.: the above expression forNOTEρis simplified because we are designing zero-mean signals. If the means are non-zero, they must be subtracted as shown here: [link].

Now, although it is not always easy to obtain partially correlated signals, there are many ways of producinguncorrelatedsignals (i.e. where the cross-correlation is zero). Therefore, let's start by considering two zero-mean, unit-power, Normally distributed complex random variables, s_{1}(t) and s_{2}(t), which are uncorrelated such that:

Or, in Matlab (GNU Octave) code, we simply write:

A useful property of Normally distributed random variables is that when we add them together, the result is also Normally distributed. Or, as more correctly described in this article [link]:

Therefore, let's choose xthe sum of two independent normally distributed random variables is normal, with its mean being the sum of the two means, and its variance being the sum of the two variances_{1}(t) to be:

since this satisfies Equation (3), because:

Then, we just have to design an x_{2}(t) which will satisfy Equations (4) and (5). Recalling (from the quote above) that linear combinations of Normally distributed random variables are also Normally distributed, let us consider constructing x_{2}(t) as follows:

Then, we need to find (complex) values for A_{1}and A_{2}which satisfy our remaining requirements (Equations 4 and 5):

An initial observation is that these conditions are independent of the complex phase of A_{2}. Therefore, we can assume A_{2}to be purely real (so that A_{2}=|A_{2}|). So, solving these simultaneous equations (Equations 12 and 13) for real-valued A_{2}therefore yields:

Finally, we rearrange Equation 13 to get:

Putting these results together, it follows that we can construct our two partially correlated signals, x_{1}(t) and x_{2}(t), as:

This gives us a neat result for 2 signals, but to generalise to N signals, we will need a more general approach. This will be described in Part II of this blog.

To round off this part, here is a simple Matlab (GNU Octave) example:

Matlab Example

Code Matlab M - [expand] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 close all; clear all; clc; % Number of samples to generate Nsamps = 100000; % Set the required second order statistics P1 = 2.3; P2 = 0.75; rho = 0.2 - 1j*0.3; % 2x2 (theoretical) covariance matrix Rxx = [P1, sqrt(P1*P2)*rho; sqrt(P1*P2)*conj(rho), P2] % Generate zero-mean, unit-power, uncorrelated Gaussian random signals s1 = randn(1, Nsamps); s2 = randn(1, Nsamps); % Create partially correlated signals x = sqrt(P1)*s1; y = sqrt(P2)*(conj(rho)*s1 + sqrt(1-abs(rho)^2)*s2); % 2x2 (practical) sample covariance matrix Rxx_samp = [x;y]*[x;y]' / Nsamps

The above code will produce an output similar to the following:

Note that we never expect the sample covariance matrix to exactly match the theoretical covariance matrix. Sample statistics are themselves random variables, having their own probability distributions.Code:Rxx = 2.3000 + 0.0000i 0.2627 - 0.3940i 0.2627 + 0.3940i 0.7500 + 0.0000i Rxx_samp = 2.3125 + 0.0000i 0.2687 - 0.3962i 0.2687 + 0.3962i 0.7525 - 0.0000i

]]>

I was recently asked how to estimate the frequency of a noisy sine wave over a short capture. In this context, 'short' means the total capture durationrelative to the period of the sine wave. So we may have, at most, a few periods of the sine wave (or, at worst, less than one period). Here is an example (download raw data: sine_data_40k.zip):

Clealy, there is less than one period available. To illustrate the difficulty with this problem, I will first describe a popular approach for estimating frequency over 'long' captures.

FFT-based Frequency Estimation ('Long' Captures)

Whenever we are interested in the frequency content of a signal, the Fast Fourier Transform (FFT) is often an excellent tool to use. In practice, when working with real-world (finite and imperfect) data, rather than using the raw output of the FFT, it is often preferable to enhance the FFT's estimate of Power Spectral Density (PSD) using Welch's method, since this can significantly reduce noise (at the expense of reduced resolution).

Here is some Matlab code showing how to use Welch's method to estimate frequency over a 'long' capture (download: pwelch_test.zip):

Code Matlab M - [expand] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 close all; clear all; clc; % Assume we capture Nsamps samples at 1kHz sample rate Nsamps = 32768; fsamp = 1000; Tsamp = 1/fsamp; t = (0:Nsamps-1)*Tsamp; % Assume the noisy signal is exactly 123Hz fsig = 123; signal = sin(2*pi*fsig*t); noise = 1*randn(1,Nsamps); x = signal + noise; % Plot time-domain signal subplot(2,1,1); plot(t, x); ylabel('Amplitude'); xlabel('Time (secs)'); axis tight; title('Noisy Input Signal'); % Choose FFT size and calculate spectrum Nfft = 8192; [Pxx,f] = pwelch(x,gausswin(Nfft),Nfft/2,Nfft,fsamp); % Plot frequency spectrum subplot(2,1,2); plot(f,Pxx); ylabel('PSD'); xlabel('Frequency (Hz)'); grid on; % Get frequency estimate (spectral peak) [~,loc] = max(Pxx); FREQ_ESTIMATE = f(loc); title(['Frequency estimate = ',num2str(FREQ_ESTIMATE),' Hz']);

The code outputs this figure:

Although the input looks like a noisy mess, we can see an extremely narrow spike in the spectrum at 123.0469 Hz (i.e. a tiny error of 0.0469 Hz).

However, if we try the same approach with the 'short' data capture (using this code: short_pwelch.zip), we get:

Now the frequency spike is invisible because it is squeezed right down in the first bin after DC. The reason is that, even though I have eliminated any averaging in the Welch method, the frequency bin spacing (i.e. resolution) is still only:

Therefore, even though the input signal is somewhere in the region of 35kHz, it is impossible for the FFT to provide an answer anywhere between 0 and 50kHz.

It is worth reiterating that the 'short' capture is short in terms of its total duration compared to the period of the sine wave. In fact, there are significantlymoresamples in the 'short' capture than the 'long' capture, but this doesn't ultimately help.

Frequency Estimation Via Curve Fitting

A far more effective method with 'short' captures is curve fitting. In a previous blog, I wrote about polynomial curve fitting. Well, in this case, we want to do sinusoidal curve fitting. This is easy to do in Matlab, using the 'sin1' option with the fit() function (code: short_sine_fit.zip):

Code Matlab M - [expand] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 close all; clear all; clc; % Load data load sine_data_40k.dat x = sine_data_40k; % Data info Tsamp = 5e-10; fsamp = 1/Tsamp; Nsamps = length(x); t = (0:Nsamps-1).'*Tsamp; % Get sine fit F = fit(t, x, 'sin1'); % Plot data and sine fit figure(); plot(F,t,x); grid on; xlabel('Time (secs)'); ylabel('Amplitude'); title(['Estimated frequency = ', num2str(F.b1/(2*pi)), ' Hz']);

The code outputs this graph:

The frequency estimate of 35526.7 Hz is very much more accurate than the FFT-based method.

]]>

We often plot some data on a graph and wonder if there might be a simple mathematical function to describe what we see:

In many cases, apolynomial approximationwould be particularly useful. That is, we want to describe the line in the form:

for some maximum powerN(called theorderof the polynomial). There is an unbelievably simple trick for doing this, which for the above graph tells us that the data is well approximated by the following 3rd order polynomial:

We can confirm the accuracy of this approximation by plotting it with the original data:

I will now describe how to do this. I will give specific details (andsource code) for a Matlab implementation, but will also describe the theory so you can implement the method in any way you choose.

The Method

This method is so simple, it can be done using just 2 lines of Matlab code:

Code Matlab M - [expand] 1 2 C = bsxfun(@power, (0:L).', (0:N)); a = pinv(C)*y;

The first line simply creates a (L+1 x N+1) matrix, C, whose (n,p)th element is n^{p}. The second line forms the Moore-Penrose pseudoinverse of C and multiplies this with the input data vector, y.

Why does this work? Well, if we insert x = 0, 1, 2, ..., L into Equation (1), then we are saying that the data values will be approximated by:

We can rewrite this as a matrix equation:

or, defining a variable for each of the bracketed terms, simply:

So, all we have to do now is find the vectorawhich best fits this equation (allowing us to then simply plug a_{0}, a_{1}, a_{2}, ..., a_{N}back into Equation (1)). I say "best fits" because an exact solution may not exist. For L>N, this is anoverdeterminedsystem of equations, so we may have to settle for an approximate solution. The most popular approximate solution I know of is the "minimum norm" solution, which is given by:

whereC^{+}denotes the Moore-Penrose pseudoinverse ofC. And... that's all there is to it! We can check our approximation by substitutinga_{approx}back into Equation (3):

and plotting this againsty, as I did above.

Final Comments

There are two issues that should be mentioned:

- Overfitting
- Numerical accuracy

If polynomial order, N, is set too high, overfitting may occur. This means that obscure higher-order terms may appear in the solution which only act to describe noise (rather than the true underlying process). Always try lower values of N first.

Numerical rounding errors can become an issue, particularly for larger N. Again, stick to lower N where possible. I have also heard that a QR decomposition of C will often yield a slightly more numerically-accurate solution fora. In this case, the Matlab code for computinga_{approx}would be:

This solution can then be refined as follows:

(and this refinement can be repeated iteratively).

Code Matlab M - [expand] 1 2 3 r = y - C*a; err = R\(R'\(C'*r)); a = a + err;

Matlab Code

Full Matlab code:

poly_approx.zip

Simplified code with minimal comments:

Code Matlab M - [expand] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 % Maximum order of polynomial N = 3; % Generate some input data L = 24; % maximum x value (starting from 0) x = (0:L).'; y = (1/3)*x.^3 - 10*x.^2 + 75*x - 33; y = y + 10*randn(L+1,1); % add some random noise to make it less obvious % Calculate a (approx) C = bsxfun(@power, (0:L).', (0:N)); a = pinv(C)*y; % Check our approximation of y y_approx = C*a; % Plot results plot(x,y); grid on; hold on; plot(x,y_approx,'--r'); xlabel('x'); ylabel('y'); legend('Input data', 'Polynomial Approximation');

]]>

Introduction

I didn't want to buy an expensive variable power supply for my small, low-power projects, so decided to make my own. I couldn't find much inspiration online, so I thought I'd write up my design in detail here in case it helps anyone else..Circuit schematics, layout, source code, compiled hex file and a complete parts list are all provided

The main features are:

- 5V - 15V regulated DC output
- Powered by standard micro USB (phone charger) connector
- High efficiency "boost" switch-mode power supply
- Optional built-in output voltage meter with dual-digit display
- Extremely compact, cheap design

Before diving into the technical detail, here's the finished product:

And here is some additional detail of the layout. The "skywiring" next to the dual 7-segment display is a bit cluttered. However, full schematics will be provided below for clarification.

Project Constraints

USB chargers are everywhere these days, so I wanted this power supply to take its input from a micro USB connector. Since USB is limited to 100mA at 5V, this would be for low-current applications only (no more than 30mA at 15V). Apart from that, I wanted all other parts to be standard, cheap and easy to solder (through-hole).

I wanted a step-up variable supply able to achieve 5V - 15V output. With minimal modification, this project can be reconfigured as a step-down (providing 0V - 5V variable). For combined step-up step-down, it would need a little more thought.

The main physical constraint of this project was that I was determined to fit the whole thing on the smallest stripboard (Veroboard) type. These have a 0.1" pitch and comprise a grid of just 9x25 holes.

The Power Supply

To keep things simple, I based the boost switch-mode power supply around the cheap and widely-available MC34063 switching regulator IC. The full details of designing this circuit are outlined in this application note: [link]. However, a neat shortcut is to use this nifty design tool: [link].

The values are pretty flexible, so here is my schematic of what I built using bits I had lying around:

This is actually just a small part of the final circuit. So, if you don't need a voltage display, you can stop after you get this far:

Construction Tip

The only difficult (non-through-hole) soldering job in this project is the micro USB connector, which has tiny pins. Since the USB cable will be pulled in and out, the connection of the micro-USB socket to the stripboard needs to be mechanically strong. To achieve this, I decided to solder five standard 0.1" pitch header pins directly to the micro-USB socket casing. Header pins are convenient for doing this, since they are held in the correct alignment by the black plastic holder:

After soldering, you can slide the plastic part off, giving something like this:

Although soldering the electrical connections is still fiddly, the mechanical connections are now through-hole and very strong.

The Voltage Meter

Hardware

The voltage meter can be built around just about any microcontroller with:

- An Analogue-to-Digital Converter (ADC)
- At least 9 digital output pins for driving the dual 7-segment digit displays

I only chose to use the PIC16F88 because I managed to pick up a batch of 64 on eBay for £20 a few months ago.

In total, the microcontroller needs to drive 15 LED segments (7 + 7 + one decimal point). This can be done using only 9 output pins by harnessing the magic of multiplexing (switching between the digits so quickly that it looks like both digits are always illuminated). A good discussion of multiplexing 7-segment displays can be found here: [link].

In order to achieve more accurate voltage readings, I chose to add a simple 3.3V external voltage reference for the ADC. I did this using the super-cheap TL431 Precision Programmable Reference IC.

So, putting all that together, here's the circuit schematic:

You can use just about any dual-digit 7-segment display for this project as long as a decimal point is available. I chose the OSL20561-IX because that's what my usual supplier had in stock and it fits my stripboard perfectly. This is a "common-anode" type device. Of course, a "common-cathode" device is also fine - you just have to wire the control circuitry low-side rather than high-side.

Note: The "Digit Select" circuitry in the schematic corresponds to the slightly cluttered "skywiring" that can be seen in the photos above. The way the corresponding segments of the two digits are wired together is not visible in the photos, so here is a picture for clarification:

The letters A-F in the picture denote the segment name (see OSL20561-IX datasheet) and dp denotes decimal point.

To pack the transistors into a tight space, I placed their flat edges face-to-face (and top-to-bottom) and wrapped the connecting legs around so they hugged each other. After soldering, they then stay firmly together.

Software

There are a couple of tricks worth noting regarding the voltmeter software I wrote to run on the microcontroller.

Firstly, I used a look-up table (LUT) to make it easy to activate the correct LEDs in order to display the digits 0-9. This table summarises which segments (A-G) need to be activated to display the required digit N (0-9):

In the last column, "Byte Code" is the 8-bit number we need to store in the LUT, assuming "g" is bit 0, "f" is bit 1, "e" is bit 2 and so on. Bit 7 is reserved for the decimal point.Code:N A B C D E F G Byte Code

0 1 1 1 1 1 1 0 126 1 0 1 1 0 0 0 0 48 2 1 1 0 1 1 0 1 109 3 1 1 1 1 0 0 1 121 4 0 1 1 0 0 1 1 51 5 1 0 1 1 0 1 1 91 6 0 0 1 1 1 1 1 31 7 1 1 1 0 0 0 0 112 8 1 1 1 1 1 1 1 127 9 1 1 1 0 0 1 1 115

A second point to note is that I used basic IIR filtering to smooth the ADC input readings. This helped to prevent the voltage display from flickering between different values. I also used the (filtered) input variance to help determine when the user is actually changing the voltage intentionally. I experimented with a few approaches and stuck with values that worked for me.

I compiled my code using the free HI-TECH PICC (v9.83) compiler. The source code and compiled hex file are included in this attachment: voltmeter.zip.

Parts List

Almost all of these parts can be changed for similar equivalents. I used whatever I had or could obtain easily. My only recommendation would be to use a good quality inductor, capable of handling several hundred mA maximum current.

Power Supply

- 1x MC34063 switching regulator IC
- 1x good quality radial leaded inductor ≥ 500μH (approx.)
- Capacitors: 1x 470pF (approx.) ceramic, 1x ≥ 470μF (approx.) electrolytic
- Resistors: 1x 220Ω, 1x 1kΩ
- 1x 10kΩ potentiometer
- 1x Diode
- 1x Micro-USB connector socket
- 1x small slide switch (SPDT)
- 7x 0.1" pitch header pins

(Optional) Voltage Meter

- Any suitable microcontroller (e.g. PIC16F88)
- DIP18 socket for microcontroller
- 1x TL431 precision programmable reference IC
- Any suitable dual-digit 7-segment display (e.g. OSL20561-IX)
- 1x NPN transistor
- 1x PNP transistor
- 1x 100nF ceramic capacitor
- Resistors: 1x 6.8kΩ, 1x 4.7kΩ, 1x 2.2kΩ, 3x 1kΩ, 1x 470Ω, 8x 330Ω (approx.)
- 1x Diode
- (Optional) 2x 0.1" pitch header pins and jumper for enabling/disabling voltage display