-
Notifications
You must be signed in to change notification settings - Fork 6
/
fourierdemo.py
176 lines (148 loc) · 5.34 KB
/
fourierdemo.py
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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Fourier transform using Python (version 3.5).
This code accompanies the book: The Fourier Transform: A Tutorial Introduction by James V Stone.
https://jim-stone.staff.shef.ac.uk/Fourier
To run code, set testing to 0 to analyse data in supplied sound file or 1 to test code on synthetic data containing just two frequencies.
Note that this code is written for transparency, not speed. It does not use complex numbers. Instead the Fourier coefficients are found using sine and cosine functions. The corresponding MatLab code generated Figures 1.8a and 1.8b in the book (this Python code provides similar figures).
Spyder version: 4.1.5 None
Python version: 3.8.1 64-bit
Tested on mac OS 10.15.7 (Catalina)
@author: JimStone
Copyright: All except the sound files are public domain.
The sound files are sebinpubaudiIMG_1563.wav
# and sebinpubaudiIMG_1563.m4a.
#############################
# Simple Fourier transform key variables
#
# x signal of N sampled values
# T length of signal in seconds
# R sampling rate in samples per second
# N number of samples in signal
# nvec vector of N ample numbers
# tvec vector of N time values in seconds
# fundamental vector of N elements
# containing values from 0 to 2*pi =
# fundamental frequency
# fmax maximum frequency in Hz, chosen by user
# dt interval in seconds between samples = 1/R
# fmin lowest frequency in Hz = 1/T
# nfrequencies number of frequencies in transform
# Cs vector of nfrequencies cosine coefficients
# Ds vector of nfrequencies sine coefficients
# As vector of nfrequencies Fourier amplitudes
# freqs vector of frequency values in Hz
#############################
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
testing = 0
# set to 1 for synthetic data, and 0 for voice data.
fs = 10 # set fontsize for graphs.
###############################################
if testing: # make simple data to test program.
frequency = 500;
# insert this frequency Hz
T = 0.5; # length of signal in seconds
period = 1/frequency; # in seconds
R = 44100; # set sampling rate, samples/s
N = R*T; # number of samples in T seconds
tvec = np.arange(0, N, 1)/R;
# time vector in units of seconds
# make vector containg 1 frequency = 1/period.
x1 = 2 * np.pi * tvec/period;
# signal in which each interval of period
# seconds increases by 1.5*2*pi.
x2 = 1.5*x1; # make sinusoids
x = np.sin(x1);
x = x + 0.5*np.sin(x2);
else: # use data from sound file
# set name of sound file
fname = 'sebinpubaudiIMG_1563.wav';
# read in data from file.
R, x = wavfile.read(fname)
# set length of sound segment in seconds.
segmentseconds = 0.5;
N = int(R * segmentseconds);
# sample number at end of segmentseconds.
# set start and end indices in data samples.
xmin = 0;
xmax = xmin+N;
x = x[xmin:xmax];
T = N/R; # number of seconds in recording
# endif
###############################################
N = len(x); # number of samples.
# specify maximum frequency to be analysed here
fmax = 1000; # Hz
nvec = np.arange(0, N, 1);
# vector of N samples
tvec = np.arange(0, N, 1)/R;
# vector of N samples in units of seconds
# Plot first part of signal.
nforplot=500; # number of samples for graph.
plt.figure("Signal")
plt.clf()
plt.title("Signal",fontsize=fs)
plt.xlabel("Time (seconds)",fontsize=fs)
plt.ylabel("Amplitude",fontsize=fs)
plt.plot(tvec[0:nforplot],x[0:nforplot])
plt.show()
# make fundamental sinusoidal vector.
nvec = np.arange(0,N); # indices = 0->N-1
fundamental = nvec * 2 * np.pi / N;
# fundamental now spans up to 2pi.
plt.figure("Fundamental")
plt.clf()
plt.title("Fundamental + one harmonic",fontsize=fs)
plt.xlabel("Time (seconds)",fontsize=fs)
plt.ylabel("Amplitude",fontsize=fs)
plt.plot(tvec,np.sin(fundamental),label="fundamental")
plt.plot(tvec,np.sin(2*fundamental),label="first harmonic")
plt.legend()
plt.show()
dt = 1/R; # interval between samples.
T = N/R;
fmin = 1.0/T; # lowest frequency in Hz.
# number of frequencies
nfrequencies = np.ceil(fmax/fmin);
nfrequencies = int(nfrequencies)
# make arrays to store Fourier coefficients.
Cs = np.arange(nfrequencies) * 0;
Ds = np.arange(nfrequencies) * 0;
# make array to store frequency values.
freqs = np.arange(nfrequencies) * 0;
###############################################
# Here is where the work gets done.
for n in range(1, nfrequencies):
freq = fmin*n;
freqs[n] = freq;
# make cosine wave at frequency freq.
cosinewave = np.cos(fundamental*n);
# find inner product of sound with cosine wave.
C = x.dot(cosinewave);
Cs[n-1] = C;
# make sine wave at frequency freq.
sinewave = np.sin(fundamental*n);
# find inner product of sound with sine wave.
D = x.dot(sinewave);
Ds[n-1] = D;
# work now done.
###############################################
# find amplitude of components, normalise by N.
As = ((Cs**2 + Ds**2) ** 0.5) / N
# Plot amplitude spectrum
plt.figure("Amplitude spectrum")
plt.clf()
plt.title("Amplitude spectrum",fontsize=fs)
plt.plot(freqs,As)
plt.xlabel("Frequency (Hz)",fontsize=fs)
plt.ylabel("Amplitude",fontsize=fs)
plt.show()
print('\nNumber of frequencies = ',nfrequencies)
print('\nFourier program completed.')
###############
# END OF FILE
################