用 Matlab 进行傅立叶变换
Fourier Transform with Matlab
我的问题类似于 this post，但比 this post 更笼统，而且我认为在规范化方面存在错误，无论如何使用最新版本的 Matlab (2015)。我在 CodeReview SE 上犹豫 post 这个，如果你认为它更合适，请在评论中告诉我。
我想使用 Matlab fft 验证傅立叶变换的以下代码，因为我在网上发现了相互矛盾的信息来源，包括Matlab 帮助自己，我一直无法用某些这样的 "recipes"（包括来自 MathWorks 团队的 answers，见下文）验证 Parseval 定理，尤其是那些为实际输入提取单边光谱的定理。
例如，仅在提取正频率时在网上发现的用于解释实值信号的对称频谱的倍幅似乎是错误的（帕塞瓦尔定理失败），而似乎有必要使用Matlab 中两个系数的平方根（我不知道为什么）。有些人似乎也像 Y = fft(X)/L 一样直接对 DFT 系数进行归一化，但我认为这很混乱，应该被劝阻；振幅是 defined 作为复数 DFT 系数除以信号长度的模数，系数本身不应被除以。一旦通过验证，我打算 post 将此代码作为 GitHub.
上的要点
function [frq,amp,phi] = fourier_transform( time, vals )
% FOURIER_TRANSFORM computes the Fast Fourier Transform of a given time-series.
%
% [freq,amp,phi] = fourier_transform(time,vals)
%
% Inputs:
%
% time - Nx1 column vector of equally-spaced timepoints (if not, input will be resampled).
% vals - NxM matrix of real- or complex-valued observations at previous timepoints.
%
% Outputs:
%
% frq - column vector of frequencies in Hz
% amp - corresponding matrix of amplitudes for each frequency (row) and signal (column)
% phi - corresponding unwrapped phase for each frequency (row) and signal (column)
%
% Note:
% To compute the power from the output amplitude, you need to multiply by the number of timepoints:
% power = numel(time) * amp.^2;
%
% References:
% https://en.wikipedia.org/wiki/Discrete_Fourier_transform
% make sure input time-series is uniformly sampled
assert( iscolumn(time), 'Input time should be a column vector.' );
assert( ismatrix(vals) && size(vals,1) == numel(time), 'Input values should be a matrix with same number of rows than of timepoints.' );
if std(diff(time)) > 1e-6
warning('Input time-course does not appear to be uniformly sampled. Resampling before continuing.');
[vals,time] = resample(vals,time);
end
% sampling information
nt = numel(time);
dt = time(2)-time(1);
fs = 1/dt;
df = fs/nt;
% complex spectrum coefficients
coef = fft(vals);
% real input
if isreal(vals)
% extract one-sided spectrum (spectrum is symmetric)
nfft = floor( nt/2 + 1 ); % eg 8 -> 5, and 7 -> 4
coef = coef( 1:nfft, : );
frq = (0:nfft-1)*df;
% correct amplitude values to account for previous extraction
fac = sqrt(2);
amp = fac*abs(coef)/nt;
amp(1,:) = amp(1,:)/fac; % .. except for the DC component
if mod(nt,2) == 0
amp(end,:) = amp(end,:)/fac; % .. and for the Nyquist frequency (only if nt is even)
end
% complex input
else
% shift the spectrum to center frequencies around 0
coef = fftshift( coef );
frq = fftshift( (0:nt-1)*df );
amp = abs(coef)/nt;
end
% make sure frq is a column vector and compute phases
frq = frq(:);
phi = unwrap(angle(coef));
end
示例使用
>> fs=1e3; t=transpose(0:1/fs:10); nt=numel(t); X=rand(nt,1);
>> [frq,amp,phi] = fourier_transform( t, X );
>> sum( abs(X).^2 ) - nt*sum( amp.^2 ) % Parseval's theorem
ans =
-2.7285e-11
错误示例 1
来自 Matlab 的 own example, and posted SO:
>> Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t)); % Sinusoids plus noise
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
>> sum(abs(y).^2) - NFFT*sum(abs(Y).^2) % Parseval's theorem
ans =
-220.4804
问题及解决方法：
来自行 Y = fft(y,NFFT)/L 中的规范化。
这应该是：
>> Y = fft(y,NFFT);
Ya = abs(Y)/NFFT; % correctly normalised amplitudes
sum(abs(y).^2) - NFFT*sum(Ya.^2) % Parseval's theorem
错误示例 2
来自 MathWorks 团队自己的 clarification post:
>> Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t)); % Sinusoids plus noise
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
NumUniquePts = ceil((NFFT+1)/2);
Y=Y(1:NumUniquePts);
MX=2*abs(Y);
MX(1)=MX(1)/2; % DC component
if ~rem(NFFT,2) % when NFFT is even, Y(1+Y/2) is the Nyquist frequency component of x, and needs to make sure it is unique.
MX(length(MX))=MX(length(MX))/2;
end
>> sum( abs(y).^2 ) - NFFT*sum( MX.^2 )
ans =
-5.3812e+03
问题及解决方法：
再次规范化。将 Y = fft(y,NFFT)/L; 替换为 Y = fft(y,NFFT)，并假设将 MX=2*abs(Y); 替换为 MX=2*abs(Y)/NFFT;。但是这里出现了振幅倍增的问题；校正因子似乎是 sqrt(2) 而不是 2.
错误示例 3
在 MatlabCentral 上作为 answer 找到：
>> Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
x = 0.7*sin(2*pi*Fs/8*t) + sin(2*pi*Fs/4*t);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
>> sum( abs(x).^2 ) - NFFT*sum( abs(Y).^2 )
ans =
-36.1891
问题及解决方法：
和第一个例子一样，归一化问题。改为：
Y = fft(x,NFFT);
Ya = abs(Y)/NFFT;
sum( abs(x).^2 ) - NFFT*sum( abs(Ya).^2 )
TL;DR（摘要）
很难找到 fft 使用 Matlab 正确规范化 amplitude/power 值的在线示例（例如，可以用 Parseval's theorem). This is crucial if you want to compare spectra between signals with different lengths. There is also an additional problem with real-valued signals because in that case the spectrum is often computed for positive frequency only, and therefore amplitude or power values need to be scaled to account for frequency-folding. Following the post and answers below, here is a gist 验证，我认为它可以缩放系数对于实值和复值输入正确且一致。
带回家的信息是：
- 不要直接归一化 DFT 系数（例如不要写
Y = fft(x)/L）；
- 如果您使用 n 点变换（例如
fft(x,nfft)），则归一化器是 nfft 而不是 numel(x)；
- 如果提取的是单边谱，需要根据对应的共轭对DFT系数调整amplitude/power值；
- 如果提取单边频谱，则应分别计算幅度和功率（即不要从幅度计算功率）。
振幅、功率和单边谱
如定义和解释on Wikipedia：
- DFT 系数很复杂，未 归一化，而逆 DFT 的公式在总和前带有一个
1/N 因子。这在某种意义上是自然的，因为在时间到频率方向上的移动可以看作是在具有不同频率的（正交）波的基础上的投影，而在频率到时间方向上的移动可以看作是加权叠加波浪。
- 在该叠加中，整体幅度应该是原始时间点的幅度（即它是一个反演），并且该加权平均中每个波的幅度是只需将相应的 DFT 系数的大小除以波数
|Xk| / N。同样，每波的power为|Xk|^2 / N。 Matlab 也使用该归一化（好吧，FFTW 使用）。
- 对于real-valued inputs，DFT系数是相应positive/negative频率的共轭对，除了直流分量（平均项，频率0），以及时间点数为偶数时的奈奎斯特频率。实际上，大多数应用程序通过仅针对正频率提取 DFT 系数来避免这种冗余。这导致振幅和功率的离散值变得复杂。
- 对应于成对 DFT 系数的振幅（除第一个系数和存在的奈奎斯特频率外的所有系数）可以简单地加倍，然后丢弃负频率。功率也一样。
- 与功率类似，但请注意，在这种情况下，使用调整后的幅度值计算离散功率值是 不正确的。实际上
N * amp_adjusted[k]^2 = N * (2*|Xk|/N)^2 不等于 2*|Xk|^2 / N（这是 OP 中二的平方根的来源）。因此，有必要独立地计算DFT系数的幅度和功率值（另一个不缩放它们的好理由）。
N点变换
许多在线示例使用显式 N 点变换：Y = fft(x,NFFT) 其中 NFFT 通常是 2 的幂，使 FFTW 的计算效率更高。
这种情况下的有效区别（假设NFFT >= N）是x在其末尾用0填充，直到达到NFFT时间点的长度。这意味着分解中的频率数量会发生变化，因此归一化应该相对于 NFFT 波分量进行，而不是原始的 N 时间点。
因此，网上找到的几乎所有示例在对系数进行归一化的方式上都是错误的。它不应该是 Y = fft(x,NFFT)/N，而应该是 Y = fft(x,NFFT)/NFFT —— 又一个放弃对复数系数进行归一化的习惯的好理由。
请注意，这对 Parseval 的等式没有影响，因为时域中添加的项均为零，因此它们对现在更大的总和的贡献也为零。但在频域中，添加的离散频率 将 通常对原始信号有响应，这给出了为什么获得的系数在填充和填充之间确实有很大不同的直观感受未填充的变换。
代码
因此，OP 中的代码不正确，似乎有必要同时输出振幅和功率，因为没有通用的归一化系数可以适应复杂和真实的情况，偶数或奇数的时间点。您可以找到 Gist here.
我发现这个在线 fft 示例 (https://habr.com/post/112068/) 非常有用。看看这个：
%% Parameters
Tm=5;% Length of signal (s)
Fd=512;% Sampling frequency (Hz)
Ak=0.5;% Constant component (Unit)
A1=1;% The amplitude of the first sinusoid (Unit)
A2=0.7;% Amplitude of the second sinusoid (Unit)
F1=13;% Frequency of the first sinusoid (Hz)
F2=42;% Frequency of the second sinusoid (Hz)
Phi1=0;% Initial phase of the first sinusoid (Degrees)
Phi2=37;% The initial phase of the second sinusoid (Degrees)
An=3*A1;% Noise Dispersion (Unit)
FftL=1024;% Number of Fourier Spectrum Lines
%% Generating work arrays
T=0:1/Fd:Tm;% Time Arrays
Noise=An*randn(1,length(T));% An array of random noise with a length equal to the array of time
Signal=Ak+A1*sind((F1*360).*T+Phi1)+A2*sind((F2*360).*T+Phi2);% Signal array (a mixture of 2x sinusoids and a constant component)
%% Spectral representation of the signal
FftS=abs(fft(Signal,FftL));% The amplitudes of the Fourier transform of the signal
FftS=2*FftS./FftL;% Spectrum normalization by amplitude
FftS(1)=FftS(1)/2;% The normalization of the constant component in the spectrum
FftSh=abs(fft(Signal+Noise,FftL));% The amplitudes of the Fourier transform of the signal + noise mixture
FftSh=2*FftSh./FftL;% Spectrum normalization by amplitude
FftSh(1)=FftSh(1)/2;% The normalization of the constant component in the spectrum
%% Plotting
subplot(2,1,1);
plot(T,Signal);
title('Signal');
xlabel('Time (s)');
ylabel('Amplitude (Unit)');
subplot(2,1,2);
plot(T,Signal+Noise);
title('Signal+Noise');
xlabel('Time (s)');
ylabel('Amplitude (Unit)');
F=0:Fd/FftL:Fd/2-1/FftL;% The frequency array of the calculated Fourier spectrum
figure();
subplot(2,1,1);
plot(F,FftS(1:length(F)));% Plotting of the spectrum of the Fourier signal
title('Signal spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude (Unit)');
subplot(2,1,2);
plot(F,FftSh(1:length(F)));% Plotting of the Fourier signal spectrum
title('Signal spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude (Unit)');
我的问题类似于 this post，但比 this post 更笼统，而且我认为在规范化方面存在错误，无论如何使用最新版本的 Matlab (2015)。我在 CodeReview SE 上犹豫 post 这个，如果你认为它更合适，请在评论中告诉我。
我想使用 Matlab fft 验证傅立叶变换的以下代码，因为我在网上发现了相互矛盾的信息来源，包括Matlab 帮助自己，我一直无法用某些这样的 "recipes"（包括来自 MathWorks 团队的 answers，见下文）验证 Parseval 定理，尤其是那些为实际输入提取单边光谱的定理。
例如，仅在提取正频率时在网上发现的用于解释实值信号的对称频谱的倍幅似乎是错误的（帕塞瓦尔定理失败），而似乎有必要使用Matlab 中两个系数的平方根（我不知道为什么）。有些人似乎也像 Y = fft(X)/L 一样直接对 DFT 系数进行归一化，但我认为这很混乱，应该被劝阻；振幅是 defined 作为复数 DFT 系数除以信号长度的模数，系数本身不应被除以。一旦通过验证，我打算 post 将此代码作为 GitHub.
function [frq,amp,phi] = fourier_transform( time, vals )
% FOURIER_TRANSFORM computes the Fast Fourier Transform of a given time-series.
%
% [freq,amp,phi] = fourier_transform(time,vals)
%
% Inputs:
%
% time - Nx1 column vector of equally-spaced timepoints (if not, input will be resampled).
% vals - NxM matrix of real- or complex-valued observations at previous timepoints.
%
% Outputs:
%
% frq - column vector of frequencies in Hz
% amp - corresponding matrix of amplitudes for each frequency (row) and signal (column)
% phi - corresponding unwrapped phase for each frequency (row) and signal (column)
%
% Note:
% To compute the power from the output amplitude, you need to multiply by the number of timepoints:
% power = numel(time) * amp.^2;
%
% References:
% https://en.wikipedia.org/wiki/Discrete_Fourier_transform
% make sure input time-series is uniformly sampled
assert( iscolumn(time), 'Input time should be a column vector.' );
assert( ismatrix(vals) && size(vals,1) == numel(time), 'Input values should be a matrix with same number of rows than of timepoints.' );
if std(diff(time)) > 1e-6
warning('Input time-course does not appear to be uniformly sampled. Resampling before continuing.');
[vals,time] = resample(vals,time);
end
% sampling information
nt = numel(time);
dt = time(2)-time(1);
fs = 1/dt;
df = fs/nt;
% complex spectrum coefficients
coef = fft(vals);
% real input
if isreal(vals)
% extract one-sided spectrum (spectrum is symmetric)
nfft = floor( nt/2 + 1 ); % eg 8 -> 5, and 7 -> 4
coef = coef( 1:nfft, : );
frq = (0:nfft-1)*df;
% correct amplitude values to account for previous extraction
fac = sqrt(2);
amp = fac*abs(coef)/nt;
amp(1,:) = amp(1,:)/fac; % .. except for the DC component
if mod(nt,2) == 0
amp(end,:) = amp(end,:)/fac; % .. and for the Nyquist frequency (only if nt is even)
end
% complex input
else
% shift the spectrum to center frequencies around 0
coef = fftshift( coef );
frq = fftshift( (0:nt-1)*df );
amp = abs(coef)/nt;
end
% make sure frq is a column vector and compute phases
frq = frq(:);
phi = unwrap(angle(coef));
end
示例使用
>> fs=1e3; t=transpose(0:1/fs:10); nt=numel(t); X=rand(nt,1);
>> [frq,amp,phi] = fourier_transform( t, X );
>> sum( abs(X).^2 ) - nt*sum( amp.^2 ) % Parseval's theorem
ans =
-2.7285e-11
错误示例 1
来自 Matlab 的 own example, and posted SO:
>> Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t)); % Sinusoids plus noise
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
>> sum(abs(y).^2) - NFFT*sum(abs(Y).^2) % Parseval's theorem
ans =
-220.4804
问题及解决方法：
来自行 Y = fft(y,NFFT)/L 中的规范化。
这应该是：
>> Y = fft(y,NFFT);
Ya = abs(Y)/NFFT; % correctly normalised amplitudes
sum(abs(y).^2) - NFFT*sum(Ya.^2) % Parseval's theorem
错误示例 2
来自 MathWorks 团队自己的 clarification post:
>> Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t)); % Sinusoids plus noise
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
NumUniquePts = ceil((NFFT+1)/2);
Y=Y(1:NumUniquePts);
MX=2*abs(Y);
MX(1)=MX(1)/2; % DC component
if ~rem(NFFT,2) % when NFFT is even, Y(1+Y/2) is the Nyquist frequency component of x, and needs to make sure it is unique.
MX(length(MX))=MX(length(MX))/2;
end
>> sum( abs(y).^2 ) - NFFT*sum( MX.^2 )
ans =
-5.3812e+03
问题及解决方法：
再次规范化。将 Y = fft(y,NFFT)/L; 替换为 Y = fft(y,NFFT)，并假设将 MX=2*abs(Y); 替换为 MX=2*abs(Y)/NFFT;。但是这里出现了振幅倍增的问题；校正因子似乎是 sqrt(2) 而不是 2.
错误示例 3
在 MatlabCentral 上作为 answer 找到：
>> Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
x = 0.7*sin(2*pi*Fs/8*t) + sin(2*pi*Fs/4*t);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
>> sum( abs(x).^2 ) - NFFT*sum( abs(Y).^2 )
ans =
-36.1891
问题及解决方法：
和第一个例子一样，归一化问题。改为：
Y = fft(x,NFFT);
Ya = abs(Y)/NFFT;
sum( abs(x).^2 ) - NFFT*sum( abs(Ya).^2 )
TL;DR（摘要）
很难找到 fft 使用 Matlab 正确规范化 amplitude/power 值的在线示例（例如，可以用 Parseval's theorem). This is crucial if you want to compare spectra between signals with different lengths. There is also an additional problem with real-valued signals because in that case the spectrum is often computed for positive frequency only, and therefore amplitude or power values need to be scaled to account for frequency-folding. Following the post and answers below, here is a gist 验证，我认为它可以缩放系数对于实值和复值输入正确且一致。
带回家的信息是：
- 不要直接归一化 DFT 系数（例如不要写
Y = fft(x)/L）； - 如果您使用 n 点变换（例如
fft(x,nfft)），则归一化器是nfft而不是numel(x)； - 如果提取的是单边谱，需要根据对应的共轭对DFT系数调整amplitude/power值；
- 如果提取单边频谱，则应分别计算幅度和功率（即不要从幅度计算功率）。
振幅、功率和单边谱
如定义和解释on Wikipedia：
- DFT 系数很复杂，未 归一化，而逆 DFT 的公式在总和前带有一个
1/N因子。这在某种意义上是自然的，因为在时间到频率方向上的移动可以看作是在具有不同频率的（正交）波的基础上的投影，而在频率到时间方向上的移动可以看作是加权叠加波浪。 - 在该叠加中，整体幅度应该是原始时间点的幅度（即它是一个反演），并且该加权平均中每个波的幅度是只需将相应的 DFT 系数的大小除以波数
|Xk| / N。同样，每波的power为|Xk|^2 / N。 Matlab 也使用该归一化（好吧，FFTW 使用）。 - 对于real-valued inputs，DFT系数是相应positive/negative频率的共轭对，除了直流分量（平均项，频率0），以及时间点数为偶数时的奈奎斯特频率。实际上，大多数应用程序通过仅针对正频率提取 DFT 系数来避免这种冗余。这导致振幅和功率的离散值变得复杂。
- 对应于成对 DFT 系数的振幅（除第一个系数和存在的奈奎斯特频率外的所有系数）可以简单地加倍，然后丢弃负频率。功率也一样。
- 与功率类似，但请注意，在这种情况下，使用调整后的幅度值计算离散功率值是 不正确的。实际上
N * amp_adjusted[k]^2 = N * (2*|Xk|/N)^2不等于2*|Xk|^2 / N（这是 OP 中二的平方根的来源）。因此，有必要独立地计算DFT系数的幅度和功率值（另一个不缩放它们的好理由）。
N点变换
许多在线示例使用显式 N 点变换：Y = fft(x,NFFT) 其中 NFFT 通常是 2 的幂，使 FFTW 的计算效率更高。
这种情况下的有效区别（假设NFFT >= N）是x在其末尾用0填充，直到达到NFFT时间点的长度。这意味着分解中的频率数量会发生变化，因此归一化应该相对于 NFFT 波分量进行，而不是原始的 N 时间点。
因此，网上找到的几乎所有示例在对系数进行归一化的方式上都是错误的。它不应该是 Y = fft(x,NFFT)/N，而应该是 Y = fft(x,NFFT)/NFFT —— 又一个放弃对复数系数进行归一化的习惯的好理由。
请注意，这对 Parseval 的等式没有影响，因为时域中添加的项均为零，因此它们对现在更大的总和的贡献也为零。但在频域中，添加的离散频率 将 通常对原始信号有响应，这给出了为什么获得的系数在填充和填充之间确实有很大不同的直观感受未填充的变换。
代码
因此，OP 中的代码不正确，似乎有必要同时输出振幅和功率，因为没有通用的归一化系数可以适应复杂和真实的情况，偶数或奇数的时间点。您可以找到 Gist here.
我发现这个在线 fft 示例 (https://habr.com/post/112068/) 非常有用。看看这个：
%% Parameters
Tm=5;% Length of signal (s)
Fd=512;% Sampling frequency (Hz)
Ak=0.5;% Constant component (Unit)
A1=1;% The amplitude of the first sinusoid (Unit)
A2=0.7;% Amplitude of the second sinusoid (Unit)
F1=13;% Frequency of the first sinusoid (Hz)
F2=42;% Frequency of the second sinusoid (Hz)
Phi1=0;% Initial phase of the first sinusoid (Degrees)
Phi2=37;% The initial phase of the second sinusoid (Degrees)
An=3*A1;% Noise Dispersion (Unit)
FftL=1024;% Number of Fourier Spectrum Lines
%% Generating work arrays
T=0:1/Fd:Tm;% Time Arrays
Noise=An*randn(1,length(T));% An array of random noise with a length equal to the array of time
Signal=Ak+A1*sind((F1*360).*T+Phi1)+A2*sind((F2*360).*T+Phi2);% Signal array (a mixture of 2x sinusoids and a constant component)
%% Spectral representation of the signal
FftS=abs(fft(Signal,FftL));% The amplitudes of the Fourier transform of the signal
FftS=2*FftS./FftL;% Spectrum normalization by amplitude
FftS(1)=FftS(1)/2;% The normalization of the constant component in the spectrum
FftSh=abs(fft(Signal+Noise,FftL));% The amplitudes of the Fourier transform of the signal + noise mixture
FftSh=2*FftSh./FftL;% Spectrum normalization by amplitude
FftSh(1)=FftSh(1)/2;% The normalization of the constant component in the spectrum
%% Plotting
subplot(2,1,1);
plot(T,Signal);
title('Signal');
xlabel('Time (s)');
ylabel('Amplitude (Unit)');
subplot(2,1,2);
plot(T,Signal+Noise);
title('Signal+Noise');
xlabel('Time (s)');
ylabel('Amplitude (Unit)');
F=0:Fd/FftL:Fd/2-1/FftL;% The frequency array of the calculated Fourier spectrum
figure();
subplot(2,1,1);
plot(F,FftS(1:length(F)));% Plotting of the spectrum of the Fourier signal
title('Signal spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude (Unit)');
subplot(2,1,2);
plot(F,FftSh(1:length(F)));% Plotting of the Fourier signal spectrum
title('Signal spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude (Unit)');