# Initial libraries
import pandas as pd
import scipy as spy
import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, fftfreq, fftshift
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd',
'#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']
DFS 中时域是周期性的信号，定义在整个实数域上的离散的点，但是我们只能输给计算机有限长的序列。那怎么变换到频域呢？这个时候可以把这段有限长的序列当作一个周期，做周期延拓，然后再用 DFS。所以 DFS 和 DFT 的公式没什么区别。
DFT 表达式如下：
值区间进行研究。
Matlab 中 FFT 结果并没有乘以，因此为与实际幅值对应，需要进行处理。0 频率，非零频率
下方DFT代码中已经乘以，为与实际幅值对应，F[0]保持不变，才与实际信号的均值对应，为更符合实际感官认知，当截取正半边频率时，其他频率点的幅值乘以 2 即可。
N = 20
dt = 1.0 / 10
tEnd = dt * N
time = np.arange(N) * dt
np.random.seed(32)
input = np.random.rand(N) * 0.8
input += 1 * np.sin(time * 2 * np.pi * 2.0)
print("Input Data Mean:[{:.5f}] ".format(input.mean()))
fftMycode = np.full(N, 0j)
fftMycodeAbs = np.full(N, 0.0)
fftnp = np.abs(fft(input)/N)
# Dft calculation based on defination
for k in range(N):
fft_k = 0 + 0j
for n in range(N):
fft_k += input[n] * np.exp(-2.0 * np.pi * k / N * n * 1.0j)
fftMycode[k] = fft_k/N
fftMycodeAbs[k] = abs(fft_k/N)
# 为将0Hz的幅值与mean值进行对比，对FFT结果数据进行了/N处理。
if (k <= N / 2):
print(
"DFT Result Comparison:[{:2}] numpy:{:.5f}, DFT:{:.5f}, Gain:{:.10f}"
.format(k, fftnp[k], fftMycodeAbs[k], fftMycodeAbs[k] / fftnp[k]))
Freq = np.fft.fftfreq(N, dt)
freqShift = fftshift(Freq)
fig1 = plt.figure(figsize=(12, 4))
ax1 = fig1.add_subplot(1, 2, 1, xlim=(0, tEnd))
plt.grid(ls=':')
plt.plot(time, input, '-o')
ax2 = fig1.add_subplot(1, 2, 2, xlim=(Freq.min(), Freq.max()))
plt.grid(ls=':')
plt.ylim(fftnp.min(), fftnp.max() * 1.05)
plt.plot(freqShift, fftshift(fftnp), '-s', label='NumPy')
plt.plot(freqShift, fftshift(fftMycodeAbs), '-.*', label='DFT')
ax1.set_title('Time domain Series (Fs=0.1s) N=20')
ax1.set_xlabel('Time (secs)')
ax1.set_ylabel('Sine (unit)')
ax2.set_title('DFT Result')
ax2.set_xlabel('Frequency (Hz)')
ax2.set_ylabel('FFT Amplitude (unit)')
ax2.legend()
## 2. 离散加窗后的分析过程
一次FFT分析截取1帧长度的时域信号，这1帧的长度总是有限的，因为FFT分析一次只能分析有限长度的时域信号。而实际采集的时域信号总时间很长，因此，需要将采样时间很长的时域信号截断成一帧一帧长度的数据块。这个截取过程叫做信号截断，该过程因为首尾值不能保证完全相同，继而为利用离散傅里叶变换（DFS）时而进行周期延拓时，必然带来阶跃冲击，从而引入高频信号。
[公众号--谭祥军 模态空间 什么是泄漏？](https://mp.weixin.qq.com/s?__biz=MzI5NTM0MTQwNA==&mid=2247484164&idx=1&sn=fdaf2164306a9ca4166c2aa8713cacc5&scene=21#wechat_redirect)
当截断不为整个周期时：
def hanning_window(N):
return 0.5 - 0.5 * np.cos(2 * np.pi * np.arange(N) / (N - 1))
# 完整周期截取时，无频谱泄漏
N = 40
# 不完整周期截断，信号实际频率未发生变化，频谱明显泄漏
N2 = 30
dt = 1.0 / 20
tEnd1 = dt * N
FreqResolution = 1 / dt / N
time11 = np.arange(N) * dt
time12 = np.arange(N) * dt + tEnd1
np.random.seed(32)
k0 = 3.0
input1 = 1 * np.sin(time11 * 2 * np.pi * k0)
w = hanning_window(N)
fig0 = plt.figure(figsize=(16, 9))
axT0 = fig0.add_subplot(2, 2, 1, xlim=(0, tEnd1 * 2))
plt.grid(ls=':')
plt.plot(np.append(time11, time12), np.append(input1, input1), '-o')
plt.plot(np.array([2.0, 2.0]), np.array([-2, 2]), 'r-o')
plt.ylim(-1, 1)
axT0.set_title('Time domain Series (Fs={:}s) Sine:{:}Hz NFFT={:}'.format(
dt, k0, N))
axT0.set_xlabel('Time (secs)')
axT0.set_ylabel('Sine (unit)')
axT1 = fig0.add_subplot(2, 2, 2, xlim=(0, tEnd1 * 2))
plt.grid(ls=':')
plt.plot(np.append(time11, time12), np.append(input1, input1), '-o')
plt.plot(np.array([1.5, 1.5]), np.array([-2, 2]), 'r-o')
plt.plot(np.array([3.0, 3.0]), np.array([-2, 2]), 'r-o')
plt.ylim(-1, 1)
axT1.set_title('Time domain Series (Fs={:}s) Sine:{:}Hz NFFT={:}'.format(
dt, k0, N2))
axT1.set_xlabel('Time (secs)')
axT1.set_ylabel('Sine (unit)')
fig1 = plt.figure(figsize=(16, 9))
ax1 = fig1.add_subplot(2, 2, 1, xlim=(0, tEnd1 * 2))
ax1.set_title('Full Period Section')
plt.grid(ls=':')
plt.plot(time11, input1, '-o')
plt.plot(time12, input1, '-o')
fftnp1 = np.abs(fft(input1))
ax1.set_title('Periodic extension (Fs={:}s) Sine:{:}Hz'.format(dt, k0))
ax1.set_xlabel('Time (secs)')
ax1.set_ylabel('Sine (unit)')
fftnp1 = np.abs(fft(input1))
N2 = 30
tEnd2 = dt * N2
FreqResolution2 = 1 / dt / N2
time21 = np.arange(N2) * dt
input2 = 1 * np.sin(time21 * 2 * np.pi * k0)
fftnp2 = np.abs(fft(input2))
ax2 = fig1.add_subplot(2, 2, 2, xlim=(0, tEnd1 * 2))
ax2.set_title('Periodic extension (Fs={:}s) Sine:{:}Hz'.format(dt, k0))
ax2.set_xlabel('Time (secs)')
ax2.set_ylabel('Sine (unit)')
plt.grid(ls=':')
plt.plot(time21, input2, '-o')
plt.plot(time21 + tEnd2 * 1.0, input2, '-o')
plt.plot(time21 + tEnd2 * 2.0, input2, '-o')
ax3 = fig1.add_subplot(2, 2, 3)
ax3.set_xlabel('Frequency (Hz)')
ax3.set_ylabel('FFT Amplitude (unit)')
Freq = np.fft.fftfreq(N, dt)
freqShift = fftshift(Freq)
plt.grid(ls=':')
plt.xticks(np.linspace(-10, 10, 11))
plt.xlim(-10, 10)
plt.ylim(fftnp1.min(), 1)
plt.plot(freqShift,
fftshift(fftnp1) / (N / 2),
'-o',
color='green',
label='Res={:.2f}Hz'.format(FreqResolution))
ax3.legend()
ax4 = fig1.add_subplot(2, 2, 4)
ax4.set_xlabel('Frequency (Hz)')
ax4.set_ylabel('FFT Amplitude (unit)')
Freq2 = np.fft.fftfreq(N2, dt)
freqShift2 = fftshift(Freq2)
plt.grid(ls=':')
plt.xticks(np.linspace(-10, 10, 9))
plt.xlim(-10, 10)
plt.ylim(0, 1)
plt.plot(freqShift2,
fftshift(fftnp2) / (N2 / 2),
'-o',
color='green',
label='Res={:.2f}Hz'.format(FreqResolution2))
ax4.legend()
def hanning_window(N):
return 0.5 - 0.5 * np.cos(2 * np.pi * np.arange(N) / (N - 1))
# 加窗后明显减少频谱
N = 30
dt = 1.0 / 20
tEnd = dt * N
FreqResolution = 1 / dt / N
time1 = np.arange(N) * dt
time2 = np.arange(N) * dt + tEnd
np.random.seed(32)
k = 3.0
inputO = 1 * np.sin(time1 * 2 * np.pi * k)
w = hanning_window(N)
inputW = inputO * w
fig1 = plt.figure(figsize=(12, 6))
ax1 = fig1.add_subplot(2, 2, 1, xlim=(0, tEnd * 2))
ax1.set_title('Full Period Section')
plt.grid(ls=':')
plt.plot(time1, inputO, '-o', color='blue')
plt.plot(time2, inputO, '-o', color='green')
fftnpO = np.abs(fft(inputO))
fftnpW = np.abs(fft(inputW))
ax1.set_title('Time domain Series (Fs=0.1s) N=20')
ax1.set_xlabel('Time (secs)')
ax1.set_ylabel('Sine (unit)')
ax2 = fig1.add_subplot(2, 2, 2, xlim=(0, tEnd * 2))
ax2.set_title('Window Padded')
ax2.set_xlabel('Time (secs)')
ax2.set_ylabel('Sine (unit)')
plt.grid(ls=':')
plt.plot(time1, inputO, '-o', color='blue')
plt.plot(time1, inputW, '-o', color='red')
plt.plot(time2, inputO, '-o', color='green')
plt.plot(time2, inputW, '-o', color='cyan')
ax3 = fig1.add_subplot(2, 2, 3)
ax3.set_xlabel('Frequency (Hz)')
ax3.set_ylabel('FFT Amplitude (unit)')
Freq = np.fft.fftfreq(N, dt)
freqShift = fftshift(Freq)
plt.grid(ls=':')
plt.xticks(np.linspace(-10, 10, 11))
plt.ylim(0, 1)
plt.plot(freqShift, fftshift(fftnpO) / (N / 2), '-o',
label='Res={:.2f}Hz'.format(FreqResolution))
ax3.legend()
ax4 = fig1.add_subplot(2, 2, 4)
ax4.set_xlabel('Frequency (Hz)')
ax4.set_ylabel('FFT Amplitude (unit)')
Freq = np.fft.fftfreq(N, dt)
freqShift = fftshift(Freq)
plt.grid(ls=':')
plt.xticks(np.linspace(-10, 10, 11))
plt.ylim(0, 1)
plt.plot(freqShift, fftshift(fftnpW) / (N / 2) * 2.0, '-o',
label='Res={:.2f}Hz'.format(FreqResolution))
ax4.legend()
定义信号： 设有一个离散信号序列 ，。
计算频谱： 使用DFT公式计算频谱 ：
其中， 是虚数单位。
获取频谱幅值和相位： 将 表示为幅度 和相位 。
选择窗函数： 选择适当的窗函数 ，常见的包括矩形窗、汉宁窗、汉明窗等。
信号加窗： 将原始信号 与窗函数 相乘，得到加窗后的信号 ：
这一步是为了减小信号在两端的突变，避免频谱泄漏。
应用DFT： 对加窗后的信号 进行DFT，得到频谱 。
在时域加窗后，频域的频谱幅值可能会受到窗函数的影响而缩小。为了进行修正，可以引入归一化修正系数。
频域幅值比值定义为：
其中，
具体表达式为：
该方法实际计算的是加窗后每一个频率点恢复到加窗前的值，这样实际就丧失了加窗的意义。
该过程对于窗函数在不同频率点的恢复很难得到解析的表达式，即上述公式与相关。
对于离散系统可以通过时域相乘等于频域卷积的方法进行处理，代码过程如下。
即：该过程并不是对需要所有频率进行缩放（如果这样操作就起不到抑制频谱泄漏的作用），因此恢复时通常按照最大幅值进行恢复或按照能量进行恢复。通常恢复频谱的实际幅值时是所有的频率点按照同一系数统一恢复。
def hanning_window(N):
return 0.5 - 0.5 * np.cos(2 * np.pi * np.arange(N) / (N - 1))
N = 31
Fs = 20
dt = 1.0 / Fs
tEnd = dt * N
FreqResolution = 1 / dt / N
freqs = np.fft.fftfreq(N, dt)
freqsShift = fftshift(freqs)
freqsExt = np.fft.fftfreq(2 * N - 1, dt / 2)
time1 = np.arange(N) * dt
time2 = np.arange(N) * dt + tEnd
np.random.seed(32)
k = 3
inputO = 1 * np.sin(time1 * 2 * np.pi * k)
w = hanning_window(N)
inputW = inputO * w
fig1 = plt.figure(figsize=(16, 6))
ax1 = fig1.add_subplot(3, 2, 1, xlim=(0, tEnd * 2))
ax1.set_title('Not Full Period Section')
ax1.set_xlabel('Time (secs)')
ax1.set_ylabel('Sine (unit)')
plt.grid(ls=':')
plt.plot(time1, inputO, '--o', color=colors[0])
plt.plot(time2, inputO, '--o', color=colors[1])
plt.plot(time1, inputW, '-o', color=colors[3])
plt.plot(time2, inputW, '-o', color=colors[2])
# 卷积运算前采用fftshift处理，结果看上去更直观
fftConv = np.convolve(fftshift(fft(w)) / N,
fftshift(fft(inputO)) / N,
mode='same')
if (N % 2 == 0):
fftConv = np.append(fftConv[1:N], fftConv[0])
# 也可以不用先fftshift再卷积，直接卷积（mode=‘full’）再将后半段直接加到前半段，效果是一样的
# 但也需要注意N为奇数和偶数的差异对待
fftConvAng = np.angle(fftConv)
fftConvAbs = abs(fftConv)
fftOAbs = np.abs(fft(inputO)) / N
fftWAbs = np.abs(fft(inputW)) / N
fftwAbs = np.abs(fft(w)) / N
fftWAbsShift = fftshift(fftWAbs)
fftOAbsShift = fftshift(fftOAbs)
fftwAbsShift = fftshift(fftwAbs)
fftOAngShift = fftshift(np.unwrap(np.angle(fft(inputO))))
fftwAngShift = fftshift(np.unwrap(np.angle(fft(w))))
fftConvAbsShift = (fftConvAbs)
fftConvAngShift = (fftConvAng)
axW = fig1.add_subplot(3, 2, 2)
axW.set_xlabel('Frequency (Hz)')
axW.set_ylabel('FFT Amplitude (unit)')
plt.grid(ls=':')
plt.plot(freqsShift, fftOAbsShift, '-o', label='Original')
plt.plot(freqsShift, fftwAbsShift, '-s', label='Window')
plt.plot(freqsShift, fftConvAbsShift, '-o', label='ConvFreq')
plt.xlim(-Fs / 2, Fs / 2)
axW.legend()
axCompare = fig1.add_subplot(3, 2, 3)
axCompare.set_xlabel('Frequency (Hz)')
axCompare.set_ylabel('FFT Amplitude (unit)')
plt.grid(ls=':')
plt.plot(freqsShift, fftOAbsShift, '-o', label='Original')
plt.plot(freqsShift, fftWAbsShift, '-s', label='ProductTime')
plt.plot(freqsShift, fftConvAbsShift, '-.o', label='ConvFreq')
plt.xlim(-Fs / 2, Fs / 2)
axCompare.legend()
axPhase = fig1.add_subplot(3, 2, 4)
axPhase.set_xlabel('Frequency (Hz)')
axPhase.set_ylabel('FFT Amplitude (unit)')
plt.grid(ls=':')
plt.plot(freqsShift,
fftOAngShift,
'-o',
color=colors[0],
label='Original Phase')
plt.plot(freqsShift, fftwAngShift, '-o', color=colors[1], label='Window Phase')
plt.xlim(-Fs / 2, Fs / 2)
axPhase.legend()
axRatio = fig1.add_subplot(3, 2, 5)
axRatio.set_xlabel('Frequency (Hz)')
axRatio.set_ylabel('Windowed/Orginal')
plt.grid(ls=':')
ratio = fftConvAbsShift / fftOAbsShift
plt.plot(freqsShift, ratio, '-o', color=colors[2], label='Window Ratio')
plt.xlim(-Fs / 2, Fs / 2)
axRatio.legend()
axPhaseConv = fig1.add_subplot(3, 2, 6)
axPhaseConv.set_xlabel('Frequency (Hz)')
axPhaseConv.set_ylabel('FFT Amplitude (unit)')
plt.grid(ls=':')
plt.plot(freqsShift,
fftConvAngShift,
'-o',
color=colors[2],
label='Conv Phase')
plt.xlim(-Fs / 2, Fs / 2)
axPhaseConv.legend()
理论推导法是根据式(4.1)和式(4.2)直接求能量恢复系数和幅值恢复系数的一种方法。
其中能量恢复系数，和为信号函数和窗函数的时域表达式。
其中为幅值恢复系数，分别为谐波信号的幅值、加窗后峰值频率处的幅值。
对频谱进行修正中要遵循能量相等或者幅值相等原则。
幅值相等原则：对信号进行截断后，频谱分析的幅值恢复系数要使加窗后的峰值处的幅值与加同样长度的矩形窗时峰值频率处的幅值相等，但加矩形窗的峰值频率必须是无泄露的频率成分（即保证整周期采样截断）。
能量相等原则：在时域截断后，频谱分析的能量恢复系数应该使加窗后能量与同样长度的矩形窗的能量相等。
根据式(4.2)，有两种幅值恢复系数求法：
以汉宁(Hanning)窗为例，采用理论推导法求能量恢复系数和幅值恢复系数。
窗长归一化后的时域表达式为
频域表达式为
幅值恢复系数推导：根据幅值恢复系数求法1可知；计算，根据洛必达法则，，带入式(4.2)得。
能量恢复系数推导：将hanning窗的表达式直接带入式(4.1)，将式(4.1)的x(t)设为1。
理论推到法求出的恢复系数为无误差的理论值。但是在使用理论推导法求取系数时，窗函数要满足一定的条件：第一，求能量恢复系数时，窗函数的时域表达式不能过于复杂，否则无法用解析法求积分值；第二，求幅值恢复系数时，必须能根据窗函数推导出窗函数频谱的理论解析表达式，否则就不能用这种方法求解出幅值恢复系数。
窗函数不满足理论推导法所需要的条件时，可采用数值积分法和仿真计算法求解恢复系数。
求取能量恢复系数中，窗函数时域表达式过于复杂，可使用数值积分的方法。数值积分通过数值解代替解析解。根据数值积分方法和计算方式的不同，计算结果的精度也会有不同。常用的数值积分方法有梯形积分方法和辛普森(Simpson)积分方法。以复化辛普森积分为例，计算hanning窗函数的能量恢复系数，运行结果为2.6667。数值积分方法主要在于编程实现被积函数。
数值积分的方法不能计算0频率下窗函数的幅值，所以不能计算幅值恢复系数。
复化Simpson求积公式
d = @(x)(0.5+0.5*cos(pi*x))^2;
Res = 2/simpr1(d,-1,1,100); % 能量恢复系数
function [s] = simpr1(f,a,b,n)
% f为被积函数；
% a，b分别为积分的下限和上限；
% n是子区间的个数,必须是偶数；
% s是所求积分数值
h=(b-a)/(n);
s1=0;
s2=0;
for k=0:n/2-1
x=a+h*(2*k+1);
s1=s1+f(x);
end
for k=1:(n/2-1)
x=a+h*2*k;
s2=s2+f(x);
end
s=h*(f(a)+f(b)+4*s1+2*s2)/3;
end
使用数值积分的方法计算凯泽-贝塞尔(Kaiser)窗的能量恢复系数。可参考窗函数编程。
长度为N的Kaiser窗定义为
其中，可以通过调节调整窗函数的形状。
对的值取3进行数值积分计算，得出能量补偿系数为3.4375。
% 测试凯泽-贝塞尔窗函数的数值积分法求能量补偿系数
kw = @(x) subI0(3*pi*sqrt(1-(2*x-1)^2 ))/subI0(3*pi); % 窗函数，alpha值取3
kw2 = @(x) (subI0(3*pi*sqrt(1-(2*x-1)^2 ))/subI0(3*pi))^2; % 窗函数被积分计算，alpha值取3
energy_conpensation_factor = 1/simpr1(kw2,0,1,100); % 能量补偿系数
function [out] = subI0(in)
% 计算窗函数使用
out = 1;
for i = 1:20
out = out + ((in/2)^i/factorial(i))^2;
end
end
function [s] = simpr1(f,a,b,n)
% f为被积函数；
% a，b分别为积分的下限和上限；
% n是子区间的个数,必须是偶数；
% s是所求积分数值
h=(b-a)/(n);
s1=0;
s2=0;
for k=0:n/2-1
x=a+h*(2*k+1);
s1=s1+f(x);
end
for k=1:(n/2-1)
x=a+h*2*k;
s2=s2+f(x);
end
s=h*(f(a)+f(b)+4*s1+2*s2)/3;
end
选择窗函数： 选择加窗时使用的窗函数 。
计算窗函数的归一化系数： 计算窗函数的归一化系数 ：
其中， 是窗函数的长度。
计算幅值恢复修正系数： 计算幅值恢复修正系数 ：
这个修正系数用于将加窗引入的幅值缩放效应进行补偿，以便更准确地还原原始信号。
计算能量恢复修正系数： 计算能量恢复修正系数 ：
则能量恢复系数：
为计算方便取
则hanning窗能量恢复系数
def hanning_window(N):
return 0.5 - 0.5 * np.cos(2 * np.pi * np.arange(N) / (N - 1))
N=1024
w=hanning_window(N)
x=np.random.rand(N)*1
E_w1 = np.sum(x**2)/np.sum((x*w)**2)
print(E_w1)
E_w2 = N/np.sum((w)**2)
print(E_w2)
1). 如果截断的信号仍为周期信号，则不存在泄漏，无须加窗，相当于加矩形窗。
2). 如果信号是随机信号或者未知信号，或者有多个频率分量，测试关注的是频率点而非能量大小，建议选择汉宁窗，像LMS Test.Lab中默认加的就是汉宁窗。
3). 对于校准目的，则要求幅值精确，平顶窗是个不错的选择。
4). 如果同时要求幅值精度和频率精度，可选择凯塞窗。
5). 如果检测两个频率相近、幅值不同的信号，建议用布莱克曼窗。
6). 锤击法试验力信号加力窗，响应可加指数窗。
| 类型 | ENBW | 幅值误差/dB | 最高旁瓣/dB | 旁瓣衰减/dB/每10个倍频程 | 应用场合 |
|---|---|---|---|---|---|
| 矩形窗 | 1 | -3.92(36.3%) | -13.3 | -20 | 未知的周期信号 |
| 汉宁窗 | 1.5 | -1.42(15.1%) | -31.5 | -60 | 大多数一般用途 |
| 平顶窗 | 3.8 | <0.01(0.1%) | -93.6 | 0 | 要求幅值精度，多用于校准用 |
| Kaiser-Bessel | 1.8 | -1.02(11.1%) | -66.6 | -20 | 要求相对高的幅值精度和频率分辨率 |
% matlab 频域卷积及自编卷积运算与Matlab自带conv对比
clc, clear, close all;
Fs = 20;
N = Fs;
Ts = 1 / Fs;
time = linspace(0, 1-Ts, Fs)';
freq = linspace(0, Fs-Fs/N, N);
rng(0);
y1 = sin(time*2*pi*2) + rand(N, 1) * 0;
y2 = sin(time*2*pi*4);
fft1 = fft(y1) / N;
fft2 = fft(y2) / N;
fft12 = fft(y1.*y2);
c1 = convs(fftshift(fft1), fftshift(fft2));
c2 = conv(fftshift(fft1), fftshift(fft2), 'full');
% a=[1 1 2 1 ]
% b=[1 2]
% c21=convs(a,b)
% c22=conv(a,b)
fftConv = conv(fftshift(fft1), fftshift(fft2), 'same');
fftConvBack = real(ifft(fftshift(fftConv*N)));
absFFT12 = abs(fft12);
absFFTConv = abs(fftConv);
figure();
set(gcf, 'unit', 'centimeters', 'position', [10, 3, 20, 20]);
subplot(3, 1, 1)
hold on;
plot(time, y1)
plot(time, y2)
plot(time, y2.*y1, 'r-o')
plot(time, fftConvBack, 'b-.*')
grid on
subplot(3, 1, 2)
hold on;
plot(freq, fftshift(abs(fft1)), '-o')
plot(freq, fftshift(abs(fft2)), '-.*')
grid on
subplot(3, 1, 3)
hold on;
plot(freq, fftshift(absFFT12)/N, '-o')
plot(freq, (absFFTConv), '-.*')
grid on
function c = convs(x, y)
N1 = length(x);
N2 = length(y);
for k = 1:length(x) + length(y) - 1
result = 0;
for n = 1:N1
a = 0;
b = 0;
if (n > 0 && n <= N1)
a = x(n);
end
if ((k - n + 1) >= 1 && (k - n + 1) <= N2)
b = y(k-n+1);
end
result = result + a * b;
end
c(k) = result;
end
end