FHX代码

clc; clear; close all pathname='data.txt'; % 给定试验数据文件名

fid = fopen(pathname,'r'); % 打开文件，作文输入数据之用。 fre = fgetl(fid); % 读入第一行，作为文字串。 fre = str2double(fre); % 将文字串转换为双精度的实数。

[ch,COUNT] = fscanf(fid,'%f,%f\\n',[2,inf]); % 读入2行，不定列数的矩阵 fclose(fid);

ch = ch'; % 矩阵转置，变为2列矩阵。［m,n]=size(ch); % 下面求频响函数Hxf(w) NFFT = 2048/4;

NOver = round(NFFT*0.35);

[Pff,F1] = psd(ch(:,2),NFFT,fre,hanning(NFFT),NOver);

[Pfx,F2] = csd(ch(:,2),ch(:,1),NFFT,fre,hanning(NFFT),NOver); % 频响函数Hxf Hxf = Pfx ./Pff;

% Plot the absolute H(w) plot(F1,abs(Hxf));

title('Frequency Response Function'); xlabel('Frequency (Hz)'); ylabel('H(w)'); grid on

% Plot the real part of the H(w) figure

subplot(211); plot(F1,real(Hxf)); title('FRF Real part'); xlabel('Frequency (Hz)'); ylabel('Real(H(w))'); grid on

% Plot the image part of the H(w) subplot(212);

plot(F1,imag(Hxf)); title('FRF Image part'); xlabel('Frequency (Hz)'); ylabel('Imag(H(w))'); grid on