АЧХ с использованием БПФ в MATLAB - PullRequest
Новая
       13

АЧХ с использованием БПФ в MATLAB

3 голосов
/ 22 октября 2010

Вот сценарий: используя анализатор спектра, у меня есть входные значения и выходные значения.количество выборок составляет 32000, а частота дискретизации составляет 2000 выборок / с, а синусоида на входе равна 50 hz, на входе - ток, а на выходе - давление в фунтах на квадратный дюйм.

Как рассчитать частотную характеристику по этим данным, используя MATLAB, используя функцию FFT в MATLAB.

Мне удалось создать синусоидальную волну, которая выдает величину и фазовые углы, вот код, которыйя использовал: 

%FFT Analysis to calculate the frequency response for the raw data
%The FFT allows you to efficiently estimate component frequencies in data from a discrete set of values sampled at a fixed rate

% Sampling frequency(Hz)
Fs = 2000;   

% Time vector of 16 second
t = 0:1/Fs:16-1;   

% Create a sine wave of 50 Hz.
x = sin(2*pi*t*50);                                                       

% Use next highest power of 2 greater than or equal to length(x) to calculate FFT.
nfft = pow2(nextpow2(length(x))) 

% Take fft, padding with zeros so that length(fftx) is equal to nfft 
fftx = fft(x,nfft); 

% Calculate the number of unique points
NumUniquePts = ceil((nfft+1)/2); 

% FFT is symmetric, throw away second half 
fftx = fftx(1:NumUniquePts); 

% Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x
mx = abs(fftx)/length(x); 

% Take the square of the magnitude of fft of x. 
mx = mx.^2; 

% Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
% The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2.

if rem(nfft, 2) % odd nfft excludes Nyquist point
  mx(2:end) = mx(2:end)*2;
else
  mx(2:end -1) = mx(2:end -1)*2;
end

% This is an evenly spaced frequency vector with NumUniquePts points. 
f = (0:NumUniquePts-1)*Fs/nfft; 

% Generate the plot, title and labels. 
subplot(211),plot(f,mx); 
title('Power Spectrum of a 50Hz Sine Wave'); 
xlabel('Frequency (Hz)'); 
ylabel('Power'); 

% returns the phase angles, in radians, for each element of complex array fftx
phase = unwrap(angle(fftx));
PHA = phase*180/pi;
subplot(212),plot(f,PHA),title('frequency response');
xlabel('Frequency (Hz)')
ylabel('Phase (Degrees)')
grid on

я взял частотную характеристику из фазового графика при 90 фазовый угол, это правильный способ для расчета частотной характеристики?

как сравнитьэто ответ на значения, которые получают из анализатора?это перекрестная проверка, чтобы увидеть, имеет ли логика анализатора смысл или нет.

Ответы [ 2 ]

1 голос
/ 22 октября 2010

На первый взгляд выглядит хорошо, но вам не хватает пары вещей:

  • Вы должны применить оконную функцию к данным временной области до БПФ, см., Например, http://en.wikipedia.org/wiki/Window_function для управления окнами вообще и http://en.wikipedia.org/wiki/Hann_window для наиболее часто используемой оконной функции (Hann aka Hanning).

  • вы, вероятно, хотите построить амплитуду логарифма в дБ, а не просто необработанную величину

0 голосов
/ 22 июля 2013

Вам следует рассмотреть возможность использования функции cpsd () для расчета частотного отклика.Масштабирование и нормализация для различных оконных функций выполняется для вас.

Частота отклика будет тогда 

G = cpsd (output,input) / cpsd (input,input)

, затем возьмите angle(), чтобы получить разность фаз между входом и выходом.

В вашем фрагменте кода не упоминаются наборы входных и выходных данных.

...