# Harmonic Analysis of Transfer Function Output

This example extracts closed-form solutions for the coefficients of frequencies in an output signal. The output signal results from passing an input through an analytical nonlinear transfer function.

This example uses the following Symbolic Math Toolbox™ capabilities:

### Motivation

To motivate the solution, we take a simple element from circuit theory: an ideal diode (in forward bias operation). The current, `I`, is the output that depends exponentially on the input, `V`. Diodes have found use in creating devices such as mixers and amplifiers where understanding the harmonic structure of the output can be useful in characterizing the device and meeting design specifications.

```syms Is V Vo real; I = Is*(exp(V/Vo) - 1)```
`I = $\mathrm{Is} \left({\mathrm{e}}^{V/\mathrm{Vo}}-1\right)$`

If V is a linear combination of 2 signals at frequencies LO and RF, the nonlinear transfer function will mix LO and RF to create output with content at combinatorial harmonic frequency combinations: freqs = {LO, 2LO, RF, 2RF, LO-RF, LO-2RF,...}.

The objective of this example is to determine the coefficients of `freqs` in the output.

### Define Input Signal

The input signal is a linear combination of two cosine signals.

```syms c1 c2 t LO RF real; input = c1*cos(LO*t) + c2*cos(RF*t)```
`input = ${c}_{1} \mathrm{cos}\left(\mathrm{LO} t\right)+{c}_{2} \mathrm{cos}\left(\mathrm{RF} t\right)$`

### Define Space of Harmonic Frequency Combinations

Below, `harmCombinations` is a combinatorial combination of the integer multiples of the input frequencies `LO` and `RF`. We restrict the space of interest defined by 3 harmonics each in the `LO` and `RF` directions.

```n = 3; harmCombinations = [kron((0:n)',ones(n*2+1,1)),repmat((-n:n)',n+1,1)]; freqs = harmCombinations*[LO;RF];```

The first `n` frequencies are just the negative harmonic frequencies and are therefore redundant considering that the input signal is real.

`freqs = freqs(n+1:end)`
```freqs =  $\left(\begin{array}{c}0\\ \mathrm{RF}\\ 2 \mathrm{RF}\\ 3 \mathrm{RF}\\ \mathrm{LO}-3 \mathrm{RF}\\ \mathrm{LO}-2 \mathrm{RF}\\ \mathrm{LO}-\mathrm{RF}\\ \mathrm{LO}\\ \mathrm{LO}+\mathrm{RF}\\ \mathrm{LO}+2 \mathrm{RF}\\ \mathrm{LO}+3 \mathrm{RF}\\ 2 \mathrm{LO}-3 \mathrm{RF}\\ 2 \mathrm{LO}-2 \mathrm{RF}\\ 2 \mathrm{LO}-\mathrm{RF}\\ 2 \mathrm{LO}\\ 2 \mathrm{LO}+\mathrm{RF}\\ 2 \mathrm{LO}+2 \mathrm{RF}\\ 2 \mathrm{LO}+3 \mathrm{RF}\\ 3 \mathrm{LO}-3 \mathrm{RF}\\ 3 \mathrm{LO}-2 \mathrm{RF}\\ 3 \mathrm{LO}-\mathrm{RF}\\ 3 \mathrm{LO}\\ 3 \mathrm{LO}+\mathrm{RF}\\ 3 \mathrm{LO}+2 \mathrm{RF}\\ 3 \mathrm{LO}+3 \mathrm{RF}\end{array}\right)$```

### Taylor Expansion

To cover the frequency spectrum of interest, a Taylor series of order four for `I(V)` is sufficient.

`s = taylor(I, V, 'Order', 4)`
```s =  $\frac{\mathrm{Is} {V}^{2}}{2 {\mathrm{Vo}}^{2}}+\frac{\mathrm{Is} {V}^{3}}{6 {\mathrm{Vo}}^{3}}+\frac{\mathrm{Is} V}{\mathrm{Vo}}$```

Use an input signal combination of the `LO` and `RF` frequencies and express `f` in terms of `cos(LO*t)` and `cos(RF*t)`.

```f0 = subs(s, V, input); f = expand(f0)```
```f =  ```

Rewrite `f` in terms of single powers of cosines.

`f = combine(f, 'sincos')`
```f =  $\frac{\mathrm{Is} {{c}_{1}}^{2}}{4 {\mathrm{Vo}}^{2}}+\frac{\mathrm{Is} {{c}_{2}}^{2}}{4 {\mathrm{Vo}}^{2}}+\frac{\mathrm{Is} {c}_{1} \mathrm{cos}\left(\mathrm{LO} t\right)}{\mathrm{Vo}}+\frac{\mathrm{Is} {c}_{2} \mathrm{cos}\left(\mathrm{RF} t\right)}{\mathrm{Vo}}+\frac{\mathrm{Is} {{c}_{1}}^{2} \mathrm{cos}\left(2 \mathrm{LO} t\right)}{4 {\mathrm{Vo}}^{2}}+\frac{\mathrm{Is} {{c}_{1}}^{3} \mathrm{cos}\left(\mathrm{LO} t\right)}{8 {\mathrm{Vo}}^{3}}+\frac{\mathrm{Is} {{c}_{1}}^{3} \mathrm{cos}\left(3 \mathrm{LO} t\right)}{24 {\mathrm{Vo}}^{3}}+\frac{\mathrm{Is} {{c}_{2}}^{2} \mathrm{cos}\left(2 \mathrm{RF} t\right)}{4 {\mathrm{Vo}}^{2}}+\frac{\mathrm{Is} {{c}_{2}}^{3} \mathrm{cos}\left(\mathrm{RF} t\right)}{8 {\mathrm{Vo}}^{3}}+\frac{\mathrm{Is} {{c}_{2}}^{3} \mathrm{cos}\left(3 \mathrm{RF} t\right)}{24 {\mathrm{Vo}}^{3}}+\frac{\mathrm{Is} {c}_{1} {{c}_{2}}^{2} \mathrm{cos}\left(\mathrm{LO} t\right)}{4 {\mathrm{Vo}}^{3}}+\frac{\mathrm{Is} {{c}_{1}}^{2} {c}_{2} \mathrm{cos}\left(\mathrm{RF} t\right)}{4 {\mathrm{Vo}}^{3}}+\frac{\mathrm{Is} {c}_{1} {c}_{2} \mathrm{cos}\left(\mathrm{LO} t+\mathrm{RF} t\right)}{2 {\mathrm{Vo}}^{2}}+\frac{\mathrm{Is} {c}_{1} {c}_{2} \mathrm{cos}\left(\mathrm{LO} t-\mathrm{RF} t\right)}{2 {\mathrm{Vo}}^{2}}+\frac{\mathrm{Is} {c}_{1} {{c}_{2}}^{2} \mathrm{cos}\left(\mathrm{LO} t-2 \mathrm{RF} t\right)}{8 {\mathrm{Vo}}^{3}}+\frac{\mathrm{Is} {c}_{1} {{c}_{2}}^{2} \mathrm{cos}\left(\mathrm{LO} t+2 \mathrm{RF} t\right)}{8 {\mathrm{Vo}}^{3}}+\frac{\mathrm{Is} {{c}_{1}}^{2} {c}_{2} \mathrm{cos}\left(2 \mathrm{LO} t+\mathrm{RF} t\right)}{8 {\mathrm{Vo}}^{3}}+\frac{\mathrm{Is} {{c}_{1}}^{2} {c}_{2} \mathrm{cos}\left(2 \mathrm{LO} t-\mathrm{RF} t\right)}{8 {\mathrm{Vo}}^{3}}$```

### Extract and Display Coefficients

Get the non-constant i.e. non-DC harmonic frequency terms of the form `cos(freq*t)`.

```cosFreqs = cos(expand(freqs*t)); terms = collect(setdiff(cosFreqs', sym(1)));```

Extract the coefficients for all harmonic frequency terms including DC.

```newvars = sym('x', [1,numel(terms)]); [cx, newvarsx] = coeffs(subs(f,terms,newvars), newvars); tx = sym(zeros(1,numel(cx))); for k=1:numel(newvarsx) if newvarsx(k) ~= 1 tx(k) = terms(newvars == newvarsx(k)); else tx(k) = newvarsx(k); end end cx = simplify(cx);```

Display the coefficients using table, `T`. Use `cosFreqs` as a row identifier.

```cosFreqs = arrayfun(@char,cosFreqs,'UniformOutput',false); Frequencies = arrayfun(@char,freqs,'UniformOutput',false); Coefficients = num2cell(zeros(size(freqs))); T = table(Frequencies,Coefficients,'RowNames',cosFreqs);```

Assign `cx` to the appropriate rows of `T` corresponding to the cosine terms of `tx`.

```nonzeroCosFreqs = arrayfun(@char,tx,'UniformOutput',false).'; T(nonzeroCosFreqs,'Coefficients') = arrayfun(@char,cx,'UniformOutput',false).';```

Now, remove row names as they are redundant.

`T.Properties.RowNames = {};`

Observe that the expressions for the terms are symmetric in LO and RF.

`T`
```T=25×2 table Frequencies Coefficients _______________ _____________________________________________ {'0' } {'(Is*(c1^2 + c2^2))/(4*Vo^2)' } {'RF' } {'(Is*c2*(8*Vo^2 + 2*c1^2 + c2^2))/(8*Vo^3)'} {'2*RF' } {'(Is*c2^2)/(4*Vo^2)' } {'3*RF' } {'(Is*c2^3)/(24*Vo^3)' } {'LO - 3*RF' } {[ 0]} {'LO - 2*RF' } {'(Is*c1*c2^2)/(8*Vo^3)' } {'LO - RF' } {'(Is*c1*c2)/(2*Vo^2)' } {'LO' } {'(Is*c1*(8*Vo^2 + c1^2 + 2*c2^2))/(8*Vo^3)'} {'LO + RF' } {'(Is*c1*c2)/(2*Vo^2)' } {'LO + 2*RF' } {'(Is*c1*c2^2)/(8*Vo^3)' } {'LO + 3*RF' } {[ 0]} {'2*LO - 3*RF'} {[ 0]} {'2*LO - 2*RF'} {[ 0]} {'2*LO - RF' } {'(Is*c1^2*c2)/(8*Vo^3)' } {'2*LO' } {'(Is*c1^2)/(4*Vo^2)' } {'2*LO + RF' } {'(Is*c1^2*c2)/(8*Vo^3)' } ⋮ ```

### Verify Coefficients

As shown below, the output waveform is reconstructed from the coefficients and has an exact match with the output.

`simplify(f0 - (dot(tx,cx)))`
`ans = $0$`

### Plot Nonlinear Transfer

The following shows the particular nonlinear transfer function analyzed above, in the time and frequency domains, for certain values of the frequencies and voltage ratios. First, extract the data.

```sample_values = struct('c1',0.4,'c2',1,'LO',800,'RF',13600,'Vo',1,'Is',1); sample_input = subs(input,sample_values)```
```sample_input =  $\frac{2 \mathrm{cos}\left(800 t\right)}{5}+\mathrm{cos}\left(13600 t\right)$```
`sample_output = subs(f,sample_values)`
```sample_output =  $\frac{127 \mathrm{cos}\left(800 t\right)}{250}+\frac{\mathrm{cos}\left(1600 t\right)}{25}+\frac{\mathrm{cos}\left(2400 t\right)}{375}+\frac{\mathrm{cos}\left(12000 t\right)}{50}+\frac{\mathrm{cos}\left(12800 t\right)}{5}+\frac{233 \mathrm{cos}\left(13600 t\right)}{200}+\frac{\mathrm{cos}\left(14400 t\right)}{5}+\frac{\mathrm{cos}\left(15200 t\right)}{50}+\frac{\mathrm{cos}\left(26400 t\right)}{20}+\frac{\mathrm{cos}\left(27200 t\right)}{4}+\frac{\mathrm{cos}\left(28000 t\right)}{20}+\frac{\mathrm{cos}\left(40800 t\right)}{24}+\frac{29}{100}$```
```sample_freqs = zeros(size(tx)); for k=1:numel(tx) cosTerm = subs(tx(k),sample_values); freq = simplify(acos(cosTerm),'IgnoreAnalyticConstraints',true)/t; sample_freqs(k) = double(freq); end sample_heights = double(subs(cx,sample_values));```

Then, use `fplot` and `stem` to plot the functions and their harmonic frequencies.

```subplot(2,2,1); fplot(sample_input,[0,0.01]) title Input subplot(2,2,3); stem([sample_values.LO, sample_values.RF],[sample_values.c1,sample_values.c2]); title 'Input Frequencies' subplot(2,2,2); fplot(sample_output,[0,0.01]) title Output subplot(2,2,4); stem(sample_freqs,sample_heights) title 'Output Frequencies'```

#### Mathematical Modeling with Symbolic Math Toolbox

Get examples and videos