Problems

Problem 9.2 (PUG)

Consider the function \[\begin{aligned} f(t) = 8 \cos(t) + 6 \sin(2 t) + \sqrt{5} \cos(4 t) + 2 \sin(4 t) + \cos(6 t - \pi/2). \end{aligned}\] (a) Find the (harmonic) magnitude and (harmonic) phase of its Fourier series components. (b) Sketch its magnitude and phase spectra. Hint: no Fourier integrals are necessary to solve this problem.

Problem 9.3 (PONYO)

Consider the function with \(a>0\) \[\begin{aligned} f(t) = e^{-a|t|}. \end{aligned}\] From the transform definition, derive the Fourier transform \(F(\omega)\) of \(f(t)\). Simplify the result such that it is clear the expression is real (no imaginary component).

clear; close all

syms t a w
assume(a > 0);
assume(t,'real');
assume(w,'real');

f = exp(-a*abs(t));

F = int(f*exp(-i*w*t),t,-Inf,Inf)

Problem 9.4 (SEESAW)

Consider the periodic function \(f:\mathbb{R}\rightarrow\mathbb{R}\) with period \(T\) defined for one period as \[\begin{aligned} f(t) = a t \quad\text{for } t \in (-T/2,T/2] \end{aligned}\] where \(a, T \in \mathbb{R}\). Perform a fourier series analysis on \(f\). Letting \(a = 5\) and \(T = 1\), plot \(f\) along with the partial sum of the fourier series synthesis, the first \(50\) nonzero components, over \(t \in [-T,T]\).

syms a t T w_n n
assume(a,'real')
assume(t,'real')
assume(T>0)
assume(w_n,'real')
assume(n,{'real','integer'})

w_n = 2*pi*n/T;  % natural frequency
fT = a*t;       % function over integration interval

c_n = simplify(...
    1/T*int(fT*exp(-i*w_n*t),t,-T/2,T/2)...
)

syms N w_N
assume(N,{'positive','integer'})
w_N = 2*pi*N/T;

a_0 = 2/T*int(fT,t,-T/2,T/2)
a_n = 2/T*int(fT*cos(w_N*t),t,-T/2,T/2)
b_n = simplify(...
    2/T*int(fT*sin(w_N*t),t,-T/2,T/2)...
)

f_n = @(n_) subs(c_n*exp(i*w_n*t),n,n_);

f_N = @(N_) subs(a_n*cos(w_N*t)+b_n*sin(w_N*t),N,N_);

nmax = 50;
Sf_n = 0;
for k = 1:nmax
    Sf_n = Sf_n + f_n(k) + f_n(-k);
end

Sf_N = 0;
for k = 1:nmax
    Sf_N = Sf_N + f_N(k);
end

params.a = 5;
params.T = 1;
Sf_n_ = matlabFunction(subs(Sf_n,params));
Sf_N_ = matlabFunction(subs(Sf_N,params));

T_ = subs(T,params); % T_ is a numerical T
f = @(t_) subs(subs(fT,params),t,mod(t_-T_/2,T_)-T_/2); % mod magic
f_ = matlabFunction(f(t));

ta = linspace(-T_,T_,301);
Sf_n_a = Sf_n_(ta); % compute once (re-using to check imag components)
figure;
plot(ta,f_(ta),'b','linewidth',1.5); hold on
plot(ta,Sf_n_a,'linewidth',1.5)
xlabel('time (s)')
ylabel('f(t)')

max(imag(Sf_n_a))

figure
plot(ta,f_(ta),'b','linewidth',1.5); hold on
plot(ta,Sf_N_(ta),'linewidth',1.5)
xlabel('time (s)')
ylabel('f(t)')

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from IPython.display import display, Markdown, Latex

sp.var('t omega_n a T n',real=True);

sym_dict = {omega_n: 2*sp.pi*n/T}
cn = sp.integrate(a*t*sp.exp(-1j*omega_n*t),(t,-T/2,T/2))
display(cn.rewrite(sp.cos).simplify())

cn_num = sp.lambdify((n,a,T),cn.subs(sym_dict))

a_num = 1
T_num = 3
n_terms = 10
n_array = np.concatenate(
  (np.arange(start=-n_terms,stop=0,step=1),
   np.arange(start=1,stop=n_terms+1,step=1))
)
cn_array = cn_num(n_array,a_num,T_num)
Mn_array = abs(cn_array)
Pn_array = np.angle(cn_array)

fig, ax = plt.subplots()
plt.stem(n_array,Mn_array,use_line_collection=True)
plt.stem([0],[0],use_line_collection=True) # add n=0 harmonic
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
plt.xlabel('harmonic $n$')
plt.ylabel('harmonic amplitude');

fig, ax = plt.subplots()
plt.stem(n_array,Pn_array,use_line_collection=True)
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
plt.xlabel('harmonic $n$')
plt.ylabel('harmonic phase (rad)');

n_sum = 50
a_num = 5
T_num = 1
t_array = np.linspace(-T_num,T_num,500)
n_array = np.concatenate(
  (np.arange(start=-n_sum,stop=0,step=1),
   np.arange(start=1,stop=n_sum+1,step=1))
)
cn_array = cn_num(n_array,a_num,T_num)

en_sym = sp.exp(1j*omega_n*t).subs(sym_dict) # symbolic version
en_num = sp.lambdify((t,n,a,T),en_sym) # numerical version
# meshgrid so lambdified function can handle computing t and n
t_array_m, n_array_m = np.meshgrid(t_array,n_array)
en_array = en_num(t_array_m,n_array_m,a_num,T_num) # evaluate

f_sum = en_array.T.dot(cn_array) # transpose to match dimensions

fig, ax = plt.subplots()
plt.plot(t_array,f_sum)
plt.xlabel('time (s)')
plt.ylabel(f'partial sum of first ${n_sum}$ terms');

Problem 9.5 (TOTORO)

Consider a periodic function \(y(t)\) with some period \(T\in\mathbb{R}\) and some parameter \(A\in\mathbb{R}\) for which one period is shown in figure 9.11.

1. Perform a trigonometric Fourier series analysis of \(y(t)\) and write the Fourier series \(Y(\omega)\).

2. Plot the harmonic amplitude spectrum of \(Y(\omega)\) for \({A = T = 1}\). Consider using computing software.

3. Plot the phase spectrum of \(Y(\omega)\) for \({A = T = 1}\). Consider using computing software.

syms A n wn T t 'real' % symbolic, real

a_n1 = 4*A*int(t*cos(wn*t),t,-T/2,-T/4) + ...
2*A/T*int(cos(wn*t),t,-T/4,T/4) + ...
-4*A*int(t*cos(wn*t),t,T/4,T/2);
a_n2 = simplify(a_n1)

a_n2 =
 
(4*A*sin((T*wn)/4))/(T*wn) - 8*A*((2*cos((T*wn)/4)^2 + T*wn*cos((T*wn)/4)*sin((T*wn)/4) - 1)/wn^2 - (cos((T*wn)/4) + (T*wn*sin((T*wn)/4))/4)/wn^2)
 

p.wn = 2*pi*n/T;
subs(subs(a_n2,p),n,5)

ans =
 
(2*A)/(5*pi) + 8*A*(T^2/(40*pi) + T^2/(100*pi^2))
 

subs(subs(a_n2,p),n,0)

Error using symengine
Division by zero.

Error in sym/subs>mupadsubs (line 160)
G = mupadmex('symobj::fullsubs',F.s,X2,Y2);

Error in sym/subs (line 145)
    G = mupadsubs(F,X,Y);

a_0 = simplify(limit(subs(a_n2,p),n,0)) % ugly

a_0 =
 
-(A*(3*T^2 - 4))/4
 

v.A = 1; % given
v.T = 1; % given
C_n = matlabFunction( ...
  abs(subs( ... % for A, T
    subs(a_n2,p), ... % for wn
  v)) ...
);
C_0 = abs( ...
  subs( ... % for A, T
    subs(a_0,p), ... % for wn
  v) ...
);

n_a = 1:20;
figure
stem([0,n_a],[C_0,C_n(n_a)])
xlabel('harmonic n')
ylabel('harmonic amplitude C_n')

phase_n = zeros(size(n_a));

figure
stem([0,n_a],[0,phase_n])
xlabel('harmonic n')
ylabel('harmonic phase')

Problem 9.6 (MALL)

Consider the function \(f:\mathbb{R}\rightarrow\mathbb{R}\) defined as \[\begin{aligned} f(t) = \begin{cases} a - a |t|/T & \text{for } t \in [-T,T] \\ 0 & \text{otherwise} \end{cases} \end{aligned}\] where \(a, T \in \mathbb{R}\). Perform a fourier transform analysis on \(f\), resulting in \(F(\omega)\). Plot \(F\) for various \(a\) and \(T\).

Matlab solution

Find the Fourier transform

First, define symbolic variables and our assumptions about them.

syms a t T f(t) w
assume(a,'real')
assume(t,'real')
assume(T,{'positive','real'})
assume(w,'real')

Define the function we will transform of, \(f(t)\), as a piecewise function as

f(t) = piecewise(...
    -T < t < T, a - a*abs(t)/T, ...
    0 ...
)

\(f(t) =\left\lbrace \begin{array}{cl} a-\frac{a\,\left|t\right|}{T} & \;\textrm{if}\;\;t<T\wedge 0<T+t\\ 0 & \;\textrm{otherwise} \end{array}\right.\)

Apply the Fourier transform definition:

F = simplify(...
    int(f(t)*exp(-i*w*t),t,-Inf,Inf)...
)

\(F =-\frac{a\,{\mathrm{e}}^{-T\,w\,\mathrm{i}} \,^{T\,w\,\mathrm{i}} -1\right)}}^2 }{T\,w^2 }\)

Plotting \(F(\omega)\)

Generate cell array FF of numerically evaluatable functions with

a_ = [1,3,6];  % array of a
T_ = [1,2,3];  % array of T
for ia = 1:length(a_)
    for iT = 1:length(T_)
        FF{ia,iT} = matlabFunction(subs(F,[a,T],[a_(ia),T_(iT)]));
    end
end

Now we are ready to plot them on the same axes with

w_ = linspace(-10,10,301);
figure
for ia = 1:length(a_)
    for iT = 1:length(T_)
        plot(...
            w_,real(FF{ia,iT}(w_)),... % turns out to be all real (not guaranteed, in general!)
            'DisplayName',...          % for legend
            ['a = ',num2str(a_(ia)), ...
            ', T = ',num2str(T_(iT))], ...
            'LineWidth',1 ...
        );hold on
    end
end
legend()
hold off
xlabel('angular frequency \omega')
ylabel('F(\omega)')

The plot is shown in .

Python solution

First, load some Python packages.

from sympy import *
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display, Markdown, Latex

Define the function

Let’s define the symbolic function.

var('t',real=True)
var('a T',positive=True)
f = Piecewise(
  (a + a/T * t, And(t >= -T, t < 0)),
  (a - a/T * t, And(t >= 0, t <= T)),
  (0, And(t < -T, t > T)),
)

Verify the function by plotting.

f_lambda = lambdify((t,a,T),f,'numpy')
t_a = np.linspace(-7,7,101)

fig, ax = plt.subplots()
plt.plot(t_a,f_lambda(t_a,2,3))
plt.xlabel('harmonic $n$')
plt.ylabel('harmonic amplitude');

Transform

The function is clearly aperiodic, so the fourier transform is the correct analysis. Simply applying the fourier transform definition, we obtain the following.

var('w',real=True)
F = integrate(
  f*exp(-I*w*t),(t,-T,T)
).simplify().rewrite(cos).trigsimp().simplify()
F_lambda = lambdify((w,a,T),F,'numpy')
display(F)

\(\displaystyle \begin{cases} - \frac{2 a \left(\cos{\left(T w \right)} - 1\right)}{T w^{2}} & \text{for}\: T^{2} w^{3} \neq 0 \\T a & \text{otherwise} \end{cases}\)

Nice! Note the use of rewrite to go from exponential to trigonometric form.

Plot

Plot with matplotlib.

a_a = np.array([1,2,3])
T_a = np.array([1,2,3])
w_a = np.linspace(-5,5,200)
fig, ax = plt.subplots()
for i,a_i in enumerate(a_a):
    for j,T_j in enumerate(T_a):
        plt.plot(
          w_a,F_lambda(w_a,a_i,T_j),
          label=f'a = {a_i}, T = {T_j}'
        )
plt.legend()
plt.grid()
plt.xlabel('$\omega$ (rad/s)')
plt.ylabel('$F(\omega)$');

Problem 9.7 (MIYAZAKI)

Consider the function \(f:\mathbb{R}\rightarrow\mathbb{R}\) defined as \[\begin{aligned} f(t) = a e^{-b \abs{t - T}} \end{aligned}\] where \(a, b, T \in \mathbb{R}\). Perform a fourier transform analysis on \(f\), resulting in \(F(\omega)\). Plot \(F\) for various \(a\), \(b\), and \(T\).

Problem 9.8 (HAKU)

Consider the function \(f:\mathbb{R}\rightarrow\mathbb{R}\) defined as \[\begin{aligned} f(t) = a \cos\omega_0 t + b \sin\omega_0 t \end{aligned}\] where \(a, b, \omega_0 \in \mathbb{R}\) constants. Perform a fourier transform analysis on \(f\), resulting in \(F(\omega)\).¹

Problem 9.9 (SECRETS)

This exercise encodes a “secret word” into a sampled waveform for decoding via a discrete fourier transform (DFT). The nominal goal of the exercise is to decode the secret word. Along the way, plotting and interpreting the DFT will be important.

First, load relevant packages.

import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display, Markdown, Latex

We define two functions: letter_to_number to convert a letter into an integer index of the alphabet (a becomes 1, b becomes 2, etc.) and string_to_number_list to convert a string to a list of ints, as shown in the example at the end.

def letter_to_number(letter):
    return ord(letter) - 96

def string_to_number_list(string):
    out = [] # list
    for i in range(0,len(string)):
        out.append(letter_to_number(string[i]))
    return out # list

print(f"aces = { string_to_number_list('aces') }")

aces = [1, 3, 5, 19]

Now, we encode a code string code into a signal by beginning with “white noise,” which is broadband (appears throughout the spectrum) and adding to it sin functions with amplitudes corresponding to the letter assignments of the code and harmonic corresponding to the position of the letter in the string. For instance, the string 'bad' would be represented by noise plus the signal \[\begin{aligned} 2 \sin{2\pi t} + 1 \sin{4\pi t} + 4 \sin{6\pi t}. \end{aligned}\]

Let’s set this up for secret word 'chupcabra'.

N = 2000
Tm = 30
T = float(Tm)/float(N)
fs = 1/T
x = np.linspace(0, Tm, N)
noise = 4*np.random.normal(0, 1, N)
code = 'chupcabra' # the secret word
code_number_array = np.array(string_to_number_list(code)) 
y = np.array(noise)
for i in range(0,len(code)):
    y = y + code_number_array[i]*np.sin(2.*np.pi*(i+1.)*x)

For proper decoding, later, it is important to know the fundamental frequency of the generated data.

print(f"fundamental frequency = {fs} Hz")

fundamental frequency = 66.66666666666667 Hz

Now, we plot.

fig, ax = plt.subplots()
plt.plot(x,y)
plt.xlim([0,Tm/4])
plt.xlabel('time (s)')
plt.ylabel('$y_n$')
plt.show()

Finally, we can save our data to a numpy file secrets.npy to distribute our message.

np.save('secrets',y)

Now, I have done this (for a different secret word!) and saved the data; download it here: ricopic.one/mathematical_foundations/source/secrets.npy

In order to load the .npy file into Python, we can use the following command.

secret_array = np.load('secrets.npy')

Your job is to (a) perform a DFT, (b) plot the spectrum, and (c) decode the message! Here are a few hints.

1. Use from scipy import fft to do the DFT.

2. Use a hanning window to minimize the end-effects. See numpy.hanning for instance. The fft call might then look like

      2*fft(np.hanning(N)*secret_array)/N

   where N = len(secret_array).

3. Use only the positive spectrum; you can lop off the negative side and double the positive side.

Load some Python packages.

We define two functions: number_to_letter to convert a number (int) into a string corresponding to the index of the alphabet (1\(\rightarrow\)a, 2\(\rightarrow\)a, etc.) and number_list_to_string to convert a list of ints to a string, as shown in the examples at the end.

def number_to_letter(number):
    return chr(number+96)

def number_list_to_string(number_list):
    out = ''
    for i in range(0,len(number_list)):
        out = out + number_to_letter(number_list[i])
    return out # list

print(number_to_letter(26))
print(number_list_to_string([1, 3, 5, 19]))

z
aces

y = np.load('secrets.npy')
N = len(y)
fs = 66.66666666666667 # must be given

Perform the FFT with a Hanning window (normalization by a factor of two is due to the “window gain” that effectively divides the spectrum by two).

Y = 2*fft(np.hanning(N)*y)/N # FFT with proper normalization

We only care about the positive spectrum, so we multiply by a factor of two and take the first “half” of the spectrum, as usual.

Y_positive_zero = 2*abs(Y[range(int(N/2))])

Similarly, the frequencies can be assigned given the sample rate fs.

freq_positive_zero = np.arange(0,N/2)*fs/N

Before we try to quantify, let’s examine the spectrum, shown in , decide the length of our word, and manually read the first letter (at \(1\) Hz).

fig, ax = plt.subplots()
ax.plot(freq_positive_zero,Y_positive_zero,'r-',linewidth=2)
plt.xlabel('frequency $f$ (Hz)')
plt.ylabel('$Y_m$');
plt.show()

It looks like the length is as follows.

word_length = 11

Now to “read” each letter. First, we define a “find nearest” function (thanks: stackoverflow!).

def find_nearest_idx(array,value):
    idx = np.searchsorted(array, value, side="left")
    if (idx > 0
    and (idx == len(array) 
    or math.fabs(value - array[idx-1]) < 
    math.fabs(value - array[idx]))):
        return idx-1
    else:
        return idx

letter_number_list = []
for i in range(0,word_length):
    idx = find_nearest_idx(freq_positive_zero,i+1)
    letter_number_list.append(int(round(Y_positive_zero[idx])))

secret_word = number_list_to_string(letter_number_list)
print('the secret word is '+secret_word)

the secret word is abracadabra

Problem 9.10 (SOCIETY)

Derive a fourier transform property for expressions including function \(f:\mathbb{R}\rightarrow\mathbb{R}\) for \[\begin{aligned} f(t) \cos(\omega_0 t + \psi) \end{aligned}\] where \(\omega_0, \psi \in \mathbb{R}\).

Problem 9.11 (FLAPPER)

Consider the function \(f:\mathbb{R}\rightarrow\mathbb{R}\) defined as \[\begin{aligned} f(t) = a u_s(t) e^{-b t} \cos(\omega_0 t + \psi) \end{aligned}\] where \(a, b, \omega_0, \psi \in \mathbb{R}\) and \(u_s(t)\) is the unit step function. Perform a fourier transform analysis on \(f\), resulting in \(F(\omega)\). Plot \(F\) for various \(a\), \(b\), \(\omega_0\), \(\psi\) and \(T\).

Problem 9.12 (EASTEGG)

Consider the function \(f:\mathbb{R}\rightarrow\mathbb{R}\) defined as \[\begin{aligned} f(t) = g(t) \cos(\omega_0 t) \end{aligned}\] where \(\omega_0 \in \mathbb{R}\) and \(g:\mathbb{R}\rightarrow\mathbb{R}\) will be defined in each part below. Perform a fourier transform analysis on \(f\) for each \(g\) below for \(\omega_1 \in \mathbb{R}\) a constant and consider how things change if \(\omega_1 \rightarrow \omega_0\).

1. \(g(t) = \cos(\omega_1 t)\)

2. \(g(t) = \sin(\omega_1 t)\)

Problem 9.13 (SAVAGE)

An instrument called a “lock-in amplifier” can measure a sinusoidal signal \(A \cos(\omega_0 t + \psi) = a \cos(\omega_0 t) + b \sin(\omega_0 t)\) at a known frequency \(\omega_0\) with exceptional accuracy even in the presence of significant noise \(N(t)\). The workings of these devices can be described in two operations: first, the following operations on the signal with its noise, \(f_1(t) = a \cos(\omega_0 t) + b \sin(\omega_0 t) + N(t)\), \[\begin{aligned} f_2(t) = f_1(t) \cos(\omega_1 t) \quad\text{and}\quad f_3(t) = f_1(t) \sin(\omega_1 t). \end{aligned}\] where \(\omega_0, \omega_1, a, b \in \mathbb{R}\). Note the relation of this operation to the Fourier transform analysis of . The key is to know with some accuracty \(\omega_0\) such that the instrument can set \(\omega_1 \approx \omega_0\). The second operation on the signal is an aggressive low-pass filter. The filtered \(f_2\) and \(f_3\) are called the in-phase and quadrature components of the signal and are typically given a complex representation \[\begin{aligned} (\text{in-phase}) + j\, (\text{quadrature}). \end{aligned}\] Explain with fourier transform analyses on \(f_2\) and \(f_3\)

1. what \(F_2 = \mathcal{F}(f_2)\) looks like,

2. what \(F_3 = \mathcal{F}(f_3)\) looks like,

3. why we want \(\omega_1 \approx \omega_0\),

4. why a low-pass filter is desirable, and

5. what the time-domain signal will look like.

Problem 9.14 (STRAWMAN)

Consider again the lock-in amplifier explored in problem 9.13. Investigate the lock-in amplifier numerically with the following steps.

1. Generate a noisy sinusoidal signal at some frequency \(\omega_0\). Include enough broadband white noise that the signal is invisible in a time-domain plot.

2. Generate \(f_2\) and \(f_3\), as described in problem 9.13.

3. Apply a time-domain discrete low-pass filter to each \(f_2 \mapsto \phi_2\) and \(f_3 \mapsto \phi_3\), such as scipy’s scipy.signal.sosfiltfilt, to complete the lock-in amplifier operation. Plot the results in time and as a complex (polar) plot.

4. Perform a discrete fourier transform on each \(f_2 \mapsto F_2\) and \(f_3 \mapsto F_3\). Plot the spectra.

5. Construct a frequency domain low-pass filter \(F\) and apply it (multiply!) to each \(F_2 \mapsto F_2'\) and \(F_3 \mapsto F_3'\). Plot the filtered spectra.

6. Perform an inverse discrete fourier transform to each \(F_2' \mapsto f_2'\) and \(F_3' \mapsto f_3'\). Plot the results in time and as a complex (polar) plot.

7. Compare the two methods used, i.e. time-domain filtering versus frequency-domain filtering.