计算信号的FFT频谱并获得正确的振幅

Question

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

• 通过乘以一个因子来补偿窗口处理
• 将双边频谱转换为单边频谱 -> 将频谱幅度乘以2
• 通过将频率节点数`N`除以规范化频谱振幅

只要fft节点数等于时间信号中的节点数，它就可以正常运行。但是，当fft节点数少于时间信号节点数时，频谱振幅就会出现错误。

我建立了一个小例子，计算了三个频谱：

1. 使用带窗函数的fft()
2. 使用不带窗函数的fft()
3. 使用spectrogram()

通过更改SampleTime，可以更改时间信号的长度，通过更改N，可以更改FFT节点数。在以下设置中，所有幅度都正确，但是对于具有Hanning窗口的频谱，幅度偏差太大。为什么会这样？非常感谢您提前的任何帮助。

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()')

- Matthias La

Answer 1

你应该改变FFT的数量。

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

这里的N为1001。

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

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

将您的代码更改如下。

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])

在这段代码中，我将所有的length(y)都改为了N。并且我把ywin修改成了以下内容。

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

结果将如下所示。

- KKS