```matlab
function Hd = buttle
Fs = 1200; % 采样频率
Fpass = 0.5; % 通带频率
Fstop = 50; % 阻带频率
Apass = 1; % 通带波动(dB)
Astop = 5; % 阻带衰减(dB)
match = 'stopband'; %波段匹配
%构造滤波器对象命名方法为butter
h = fdesign.lowpass(Fpass, Fstop, Apass, Astop, Fs);
Hd = design(h, 'butter', 'MatchExactly', match);
```
一些操作,懂的都懂。
%% Set up filters.
half_rate = samp_rate / 2;
uppe_orde = 6;
uppe_stop = 40;
lowe_orde = 3;
lowe_stop = 4;
[lu,ld] = butter(uppe_orde,uppe_stop/half_rate,'low');
[hu,hd] = butter(lowe_orde,lowe_stop/half_rate,'high');
data_wind = samp_rate * 2;
iter_star = 1:samp_rate:n_data-data_wind+1;
freq_bins = [0:samp_rate/2-1 -samp_rate/2:-1]'/samp_rate;
time_bins = (1:length(iter_star) * (0.5*data_wind)) / samp_rate;
freq_sele = freq_bins <= uppe_stop / samp_rate & freq_bins >= -uppe_stop / samp_rate;
freq_bins = freq_bins(freq_sele,:);
freq_time_prof = zeros(length(freq_bins), length(time_bins),length(file_path_list));
for kk = 1:length(file_path_list),
%% Processing data streams.
chan_data = chan_data_list{kk};
chan_powe = chan_data .* conj(chan_data);
ante_pairs = zeros(2, length(iter_star));
for ii = 1:length(iter_star),
%% Select antenna pair.
part_data = chan_data(iter_star(ii):iter_star(ii)+data_wind-1,:);
part_powe = chan_powe(iter_star(ii):iter_star(ii)+data_wind-1,:);
part_vari = sqrt(var(part_powe,1));
part_mean = mean(part_powe,1);
part_vari_ante = mean(reshape(part_vari.', n_subc, n_ante));
part_mean_ante = mean(reshape(part_mean.', n_subc, n_ante));
part_indi_ante = part_mean_ante ./ part_vari_ante;
[~, maxi] = max(part_indi_ante);
[~, mini] = min(part_indi_ante);
ante_pair = [mini, maxi];
ante_pairs(:,ii) = ante_pair.';
part_diff = part_data(:,(ante_pair(1)-1)*n_subc+1:ante_pair(1)*n_subc) ...
.* conj(part_data(:,(ante_pair(2)-1)*n_subc+1:ante_pair(2)*n_subc));
% part_diff = part_data .* conj(part_data);
%% Filter out noises and interferences.
for jj = 1:n_subc,
part_diff(:,jj) = filter(lu, ld, part_diff(:,jj));
part_diff(:,jj) = filter(hu, hd, part_diff(:,jj));
end
%% PCA analysis.
pca_coef = pca(part_diff);
part_diff = part_diff * pca_coef(:,1);
part_wind = floor(0.5*samp_rate+1:1.5*samp_rate);
part_spec = tfrsp(part_diff, part_wind, samp_rate, tftb_window(samp_rate/4+1, 'gauss'));
freq_time_prof(:,(ii-1)*samp_rate+1:ii*samp_rate,kk) = part_spec(freq_sele, :);
end
freq_time_prof(:,:,kk) = abs(freq_time_prof(:,:,kk)); % ./ repmat(sum(abs(freq_time_prof(:,:,kk)),1), size(freq_time_prof(:,:,kk),1), 1);
end
% freq_time_prof = freq_time_prof(freq_bins >= 0,:,:);
% freq_bins = freq_bins(freq_bins >= 0);
save(save_path, 'freq_time_prof');
else