计算信号的 fft 频谱并获得正确的幅度

Question

我正在尝试一些相当简单的方法，即在具有固定长度的信号频谱中获取正确的值。我不知道包含频率，因此我无法相应地选择 fft 节点的数量。这给我留下了泄漏的问题。我试图通过使用窗口淡出我的时间信号的开始和结束来避免这种情况。有了这个，我在获得正确的幅度方面遇到了问题。这是我到目前为止所做的：

• 通过将频谱乘以一个因子来补偿窗口
• 双面到单面频谱 -> 将频谱幅度乘以 2
• 正态谱振幅除以频率节点数'N'

只要 fft 节点的数量等于我的时间信号中的数量，它就可以正常工作。但是当fft节点的数量少于时间信号节点时，频谱的幅度就会出错。

我构建了一个计算三个光谱的小例子：

1. 在窗口中使用 fft()
2. 在没有窗口的情况下使用 fft()
3. 使用频谱图()

通过改变 SampleTime 可以改变时间信号的长度，通过改变 N 可以改变 FFT 节点的数量。使用以下设置，所有幅度都是正确的，但对于具有汉宁窗的频谱。为什么差这么多？提前感谢您的帮助。

SampleTime=0.001;
Fs=1/SampleTime;
t=0:SampleTime:10;

y=2*sin(2*pi*10*t)+3*cos(2*pi*30*t); % Time Signal

fn = Fs/2; % Nyquistfrequency
N = 1001; % FFT-length (use N=2^x to use DFT-Algorithm)
df = Fs/N; % frequency resolution

%% Generating spectrum using fft()
win = hann(length(y))';  % generate Hanning window for Time Signal

ywin = y.*win*length(win)/sum(win); % amplitude correction to get correct values with window

Hwin = fft(ywin, N); % spectrum of windowed signal
H = fft(y,N);  % complex spectrum of signal without window

amplHwin = abs(Hwin); % absolute values of spectrum of windowed signal
amplH = abs(H); % absolute values of spectrum without window

% Double sided spectrum -> single sided spectrum of windowed signal
amplitudengang = fftshift(amplH/N);
amplitudengang = amplitudengang(ceil(length(amplitudengang)/2+1):end);
amplitudengang(2:end) = amplitudengang(2:end)*2;

% Double sided spectrum -> single sided spectrum of unwindowed signal
amplitudengangwin = fftshift(amplHwin/N);
amplitudengangwin = amplitudengangwin(ceil(length(amplitudengangwin)/2+1):end);
amplitudengangwin(2:end) = amplitudengangwin(2:end)*2;

%% Generating spectrum using spectrogram()
noverlap = floor(0.5*N);
window = hann(N);
[B,F,T] = spectrogram(y,window,noverlap,N,Fs);
% Scaling
FFTScalingFactor = 2 * 1 / N; % 2 for Hanning window, 2^(-1/2) for RMS, 1 for peak
A = FFTScalingFactor * B;
B = FFTScalingFactor * abs(B);
B(2:end-1,:) = 2.*B(2:end-1,:); % single side spectrum with even block length

%% Plots

x_fn = 0 : df : fn-df; % Frequency vector

figure
subplot(2,2,1)
plot(x_fn,amplitudengangwin,'r')
hold on
plot(x_fn,amplitudengang)
legend('Hanning window','no window')
title('FFT()')
subplot(2,2,3)
plot(t,y)
title('Time signal')
subplot(2,2,2)
plot(F,B(:,1))
title('Spectrogram()')

Answer 1

您应该更改 FFT 的数量。

Hwin = fft(ywin, N); % spectrum of windowed signal
H = fft(y,N);  % complex spectrum of signal without window

这里 N 是 1001。

这意味着当您执行窗口信号的 FFT 时，您只使用了所有数据样本中的前 1001 个样本（10001 个样本）。

您使用的窗口信号如下所示。

如下更改您的代码

Hwin = fft(ywin, length(y)); % spectrum of windowed signal
H = fft(y,length(y));  % complex spectrum of signal without window
    amplHwin = abs(Hwin); % absolute values of spectrum of windowed signal
amplH = abs(H); % absolute values of spectrum without window

% Double sided spectrum -> single sided spectrum of windowed signal
amplitudengang = fftshift(amplH/length(y));
amplitudengang = amplitudengang(ceil(length(amplitudengang)/2+1):end);
amplitudengang(2:end) = amplitudengang(2:end)*2;

% Double sided spectrum -> single sided spectrum of unwindowed signal
amplitudengangwin = fftshift(amplHwin/length(y));
amplitudengangwin = amplitudengangwin(ceil(length(amplitudengangwin)/2+1):end);
amplitudengangwin(2:end) = amplitudengangwin(2:end)*2;

fftAll=[amplitudengang' amplitudengangwin'];

figure(1)
    plot(fftAll)
    xlim([0 500])

你会得到低于频谱。

您定义了 N=1001，但您的数据长度为 10001。如果您只想使用 1001 个样本，则应将窗口长度设为 1001。

如果您只想在数据中使用 1001 个样本，请更改您的代码，如下所示。

SampleTime=0.001;
Fs=1/SampleTime;
t=0:SampleTime:10;

y=2*sin(2*pi*10*t)+3*cos(2*pi*30*t); % Time Signal

fn = Fs/2; % Nyquistfrequency
N = 1001; % FFT-length (use N=2^x to use DFT-Algorithm)
df = Fs/N; % frequency resolution

%% Generating spectrum using fft()
win = hann(N)';  % generate Hanning window for Time Signal

figure(10)
plot(win)

ywin = y(1:N).*win*length(win)/sum(win); % amplitude correction to get correct values with window
figure(11)
plot(ywin)
xlim([0 1001])
   Hwin = fft(ywin, N); % spectrum of windowed signal
H = fft(y,N);  % complex spectrum of signal without window
    amplHwin = abs(Hwin); % absolute values of spectrum of windowed signal
amplH = abs(H); % absolute values of spectrum without window

% Double sided spectrum -> single sided spectrum of windowed signal
amplitudengang = fftshift(amplH/N);
amplitudengang = amplitudengang(ceil(length(amplitudengang)/2+1):end);
amplitudengang(2:end) = amplitudengang(2:end)*2;

% Double sided spectrum -> single sided spectrum of unwindowed signal
amplitudengangwin = fftshift(amplHwin/N);
amplitudengangwin = amplitudengangwin(ceil(length(amplitudengangwin)/2+1):end);
amplitudengangwin(2:end) = amplitudengangwin(2:end)*2;

fftAll=[amplitudengang' amplitudengangwin'];

figure(1)
    plot(fftAll)
    xlim([0 50])

在这段代码中，我将所有长度（y）都更改为 N。 我改变了 ywin 如下所示。

ywin = y(1:N).*win*length(win)/sum(win);

结果将如下所示。