Réponse en fréquence utilisant FFT dans MATLAB

Question

Voici le scénario: en utilisant un analyseur de spectre, j'ai les valeurs d'entrée et les valeurs de sortie. le nombre d'échantillons est 32000 et la fréquence d'échantillonnage est 2000 échantillons/s, et l'entrée est une onde sinusoïdale de 50 hz, l'entrée est le courant et la sortie est la pression en psi. Comment calculer la réponse en fréquence à partir de ces données en utilisant MATLAB, en utilisant la fonction FFT dans MATLAB.

i a pu générer une onde sinusoïdale, qui donne les Magnitude et les angles de phase, est le code ici que je l'ai utilisé:

%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 

j'ai pris la réponse en fréquence de la courbe de phase à 90 degré angle de phase, est-ce la bonne façon de calculer la réponse en fréquence?

comment comparer cette réponse aux valeurs obtenues à partir de l'analyseur? C'est une vérification croisée pour voir si la logique de l'analyseur a un sens ou non.

Answer 1

Attend OK à première vue, mais deux ou trois choses que vous êtes absent:

• vous devez appliquer une fonction de fenêtre pour les données de domaine temporel avant la FFT, voir par exemple http://en.wikipedia.org/wiki/Window_function pour la fenêtrage en général et http://en.wikipedia.org/wiki/Hann_window pour la fonction de fenêtre la plus couramment utilisée (Hann aka Hanning).

• vous voulez probablement tracer log amplitude en dB plutôt que magnitude brute

Answer 2

Vous devriez envisager de regarder la fonction cpsd() pour le calcul de la réponse en fréquence. La mise à l'échelle et la normalisation de diverses fonctions de fenêtre sont gérées pour vous.

la reponse en fréquence serait alors

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

puis prendre la angle() pour obtenir la différence de phase entre l'entrée et la sortie.

Votre extrait de code ne mentionne pas les données d'entrée et de sortie.