Фильтрация сигналов в octave

Завершаю серию заметок по работе в octave заметкой о фильтрации сигналов. Используется такой же сигнал, который использован в предыдущих заметках, но в данном примере показана возможность использования фильтра.

set (0, 'defaulttextfontname', 'Terminus,16');

QUALITY = 20;
N = 512 ; % Максимальное значение времени t
minfreq=100/(N-1);
maxfreq = 257/(N-1); % максимальная частота волны
dfreq = 2 * maxfreq; % Частота дискретизации (замеров)
disp('Частота волны freq = '), disp(maxfreq)

t = 0:1/(QUALITY*dfreq):N;

% Уравнение волны
F = N*log(N+1)/(N^2-1)*sin((2*pi*100/N)*t+N/10) + \
(N+3)*log(N/3)/(N^2-1)*sin((2*pi*157/N)*t+N/3) + \
(N+4)*log(N/4)/(N^2-1)*cos((2*pi*257/(N-1))*t+N/4);

%%============================================================================%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NOISE (ШУМ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%============================================================================%%

figure(1)
NoiseKoef = 100;
noise=randn([1, length(F)])/NoiseKoef;
FN=F + noise; % signal with noise
subplot(3,1,1);
plot(t, F,'b', t, FN, 'r');
title('signal with noise');
axis([0, 50])
subplot(3,2,3); 

% Expected value (математическое ожидание)
ExpVal = mean(FN)

% Variance (дисперсия)
Var = var(FN)

% Correlation (корреляция)
Cor = xcorr(FN);
subplot(3,2,6)
plot(Cor,'r')
title('Correlation')
grid on

% Probability density (плотность рапределения)
subplot(3,2,3)
ProbDens = hist(FN);
hist(FN)
title('Probability density')
grid on

% Probability distribution (закон распрделения)
subplot(3,2,4)
ProbDist = 0;
for i=1:length(ProbDens)
    ProbDist(i+1)=ProbDens(i)+ProbDist(i);
end
ProbDist=ProbDist/length(FN);
plot(ProbDist);
title('Probability distribution')
grid on

% Spectrum (спектр)
subplot(3,2,5)
spectrum = abs(fft(FN));
plot(spectrum,'r')
title('Spectrum')
grid on

%%============================================================================%%
%%%%%%%%%%%%%%%%%%% Low-pass filter (фильтр низких частот) %%%%%%%%%%%%%%%%%%%%%
%%============================================================================%%
figure(2);
% зашумленный сигнал
subplot(4,1,1);
plot(t, F,'b', t, FN, 'r');
title('signal with noise');
axis([0, 50])

% исходный сигнал
subplot(4,1,2);
plot(t, F,'m');
title('original signal');
axis([0, 50])

% low-pass filter (фильтр Баттерворта)
%Wp = minfreq*0.3; Ws = maxfreq;
%[n,Wn] = buttord(Wp,Ws,3,60)
%[b,a] = butter(n,Wn);
%FLP = filter(b,a,FN);

%subplot(3,1,3)
%plot(t, FLP, 'b')
%title('low-pass filter');
%axis([0, 50])

Np=3;%	Порядок фильтра
fc=minfreq*0.3;	% cutoff frequency (частота среза)
[b,a]=butter(Np,fc);
FLP = filter(b,a,FN);

subplot(4,1,3)
plot(t, FLP, 'b')
title('low-pass filter');
axis([0, 50])

% выделение гармоники
%figure(3)
subplot(4,1,4);
spectrum = abs(fft(FLP));
plot(spectrum,'r');
title('harmonic');
axis([0, 50, 0, 5])

~ от aleos на 24 Май 2008.