Zespoły filtrów i transformacje ortogonalne
Autorzy: Marek Parfieniuk
Pliki: rozdzial_06.zip
Przykład 6.1
Przetwarzanie sygnału za pomocą funkcji MATLAB, które stanowią bezpośrednią implementację części analizującej oraz syntezującej zespołu filtrów równomiernych (plik przyklad_6_01.m).% Przykład 6.1 function przyklad_06_01() M = 8; % Liczba kanałów zespołu filtrów. h = dctIImtx(M); % Macierz współczynników filtrów analizy (wiersz odpowiada ... % filtrowi). Np. filtry określone macierzą DCT. f = fliplr(h); % Współczynniki filtrów syntezy. W przypadku DCT uzyskiwane ... % z odwrócenia kolejności współczynników filtrów analizy. S = M; % Współczynnik podpróbkowania. Np. próbkowanie krytyczne. x = [1 zeros(1, 32 - 1)]; % Sygnał wejściowy. Np. impuls jednostkowy. xs = fb_analysis(h, S, x); % Analiza sygnału wejściowego. % Sygnały kanałów w wierszach $xs$. ys = xs; % Brak przetwarzania kanałów, aby sprawdzić czy zespół filtrów zapewnia PR. y = fb_synthesis(f, S, ys); % Synteza wyjściowego sygnału sygnału pełnopasmowego. fb_show_io(x, y, 512); % Obrazowanie sygnału oryginalnego i zrekonstruowanego. fb_show_response(h, 256); % Obrazowanie charakterystyk filtrów.
Przykład 6.2
Projektowanie filtra do decymacji/interpolacji oraz demonstracja zmiany tempa próbkowania (plik przyklad_6_02.m).% Przykład 6.2:
function przyklad_06_02()
N = 32; M = 4; % Długość filtru i współczynnik decymacji/interpolacji.
f = 1 / M; % Znormalizowana częstotliwość odcięcia.
delta = 1 / (2 * M); % Znormalizowana szerokość pasma przejściowego.
fp = f - delta; fst = f; % Częstotliwości graniczne pasma przepustowego i ...
%fp = f - delta/2; fst = f + delta/2; % zaporowego: dwie alternatywne specyfikacje.
d = fdesign.lowpass('n,fp,fst', N, fp, fst); % Określenie specyfikacji filtru.
h = design(d, 'firls'); % Projektowanie filtru metodą najmniejszych kwadratów.
L = 128; % Liczba wyznaczanych punktów segmentu charakterystyki amplitudowej.
[H, w] = freqz(h, 2 * L * M, 'whole'); % Wyznaczenie charakterystyki ...
figure; hold on; % i jej prezentacja.
for k = 1 : M, % Iteracja po segmentach charakterystyki.
idx = (1 : (L + 1)) + (k - 1) * L; % Indeksy punktów segmentu.
if (mod(k, 2) == 0), idx = fliplr(idx); end; % Co drugi segment "zawijamy".
plot(w(1 : (L + 1)) / 2 / pi, 20 * log10(abs(H(idx)))); % Wyświetlenie segmentu.
end
box on; grid on; axis([0 0.5/M -100 10]);
xlabel('\Omega/2\pi'); ylabel('|H(e^{j\Omega}|) [dB]');
t = (0 : 127) / 128; % Generacja sygnału testowego ...
x = sin(t * 2 * pi * 8) + sin(t * 2 * pi * 24); % np. sumy dwóch sinusoid.
clc
if (strcmpi(input('Wybierz: D = decymacja, I = interpolacja.\n', 's'), 'D')),
t = sprintf('Decymacja: M = %d, ', M);
if (strcmpi(input(['Czy użyć filtru do eliminacji nakładania widmowego? ' ...
'T = tak, inny znak = nie.\n'], 's'), 'T')),
v = filter(h, x); t = [t 'Z filtrem']; else v = x; t = [t 'Bez filtru'];
end
y = v(1 : M : end);
else
t = sprintf('Interpolacja: M = %d, ', M);
v = zeros(1, length(x) * M); v(1 : M : end) = x;
if (strcmpi(input(['Czy użyć filtru do eliminacji replik widma? ' ...
'T = tak, inny znak = nie.\n'], 's'), 'T')),
y = filter(h, v); t = [t 'Z filtrem']; else y = v; t = [t 'Bez filtru'];
end
end
figure % Prezentacja przebiegu zmiany tempa próbkowana.
L = 128; f = (0 : (L - 1)) / L;
subplot(3, 2, 1); stem(x); ylabel('x[n]'); title(t);
subplot(3, 2, 2); stem(f, abs(fft(x, L))); ylabel('|X(e^{j\Omega})|');
subplot(3, 2, 3); stem(v); ylabel('v[n]');
subplot(3, 2, 4); stem(f, abs(fft(v, L))); ylabel('|V(e^{j\Omega})|');
subplot(3, 2, 5); stem(y); ylabel('y[n]'); xlabel('n');
subplot(3, 2, 6); stem(f, abs(fft(y, L)));
ylabel('|Y(e^{j\Omega})|'); xlabel('\Omega/2\pi');
Przykład 6.3
Generowanie i badanie filtrów modulowanych DFT (plik przyklad_6_03.m).% Przykład 6.3:
function przyklad_06_03()
M = 8; % Liczba kanałów.
S = 8; % Współczynnik podpróbkowania. Np. próbkowanie krytyczne.
N = 32; % Długość filtru prototypowego.
fc = 1/M + 0.0018; % Znormalizowania częstotliwość odcięcia filtru prototypowego.
p = fir1(N - 1, fc, kaiser(N,5)); % Synteza prototypu z wykorzystaniem okna Kaisera.
% p = ones(1, M) ./ M; % Alternatywnie: okno prostokątne jako prototyp.
hs = fb_dft_modulated_filters(p, M); % Generacja współczynników filtrów ...
% na podstawie prototypu.
fb_show_response(hs, 512); % Prezentacja charakterystyk amplitudowych.
fb_show_prototype(p, 512, 1/M); % Prezentacja prototypu.
x = [1 zeros(1, 127)]; % Sygnał wejściowy.
xs = fb_analysis(hs, S, x); % Analiza podpasmowa. Sygnały kanałów w wierszach $xs$.
ys = xs; % Brak przetwarzania kanałów, aby sprawdzić czy zespół filtrów zapewnia PR.
y = fb_synthesis(hs * S, S, ys); % Synteza wyjściowego sygnału pełnopasmowego.
fb_show_io(x, real(y), 512); % Prezentacja sygnału oryginalnego i ...
max(abs(imag(y))) % zrekonstruowanego. Część urojona wyniku powinna być zerowa.
%===================================================================================
function h = fb_dft_modulated_filters(p, M) % Współczynniki $M$-kanałowego ...
% zespołu filtrów modulowanych DFT.
N = length(p); % Długość filtru prototypowego $p$.
[k, n] = ndgrid(0 : (M - 1), 0 : (N - 1)); % Generacja indeksów i ...
h = repmat(p, M, 1) .* exp(-1i * 2 * pi / M) .^ (-k .* n); % współczynników filtrów.
Przykład 6.4
Projektowanie zespołu filtrów z modulacją kosinusową, który aproksymuje perfekcyjną rekonstrukcję (plik przyklad_6_04.m).% Przykład 6.4:
function przyklad_06_04()
rand('state', sum(100 * clock)); % Inicjalizacja generatora liczb losowych.
phi_start = 0; % Ustalenie rozwiązania początkowego.
options = optimset('MaxFunEvals', 2000, ... % Parametry algorytmu optymalizacyjnego.
'DiffMinChange', 0.00001, 'DiffMaxChange', 0.15, 'TolFun', 1e-10, 'TolX', 1e-10);
phi_opt = fminsearch(@obj_fun, phi_start, options); % Optymalizacja bez ograniczeń.
obj_fun(phi_opt); % Prezentacja najlepszego znalezionego rozwiązania.
%===================================================================================
function err = obj_fun(phi) % Funkcja celu optymalizacji.
% Zamiast bezpośrednio optymalizować częstotliwość ...
% odcięcia optymalizujemy jej odchylenie od wartości nominalnej.
% Dzięki temu nie trzeba uwzględniać ograniczeń podczas optymalizacji.
M = 32; N = 512; % Liczba kanałów i długość filtru prototypowego.
fc = (1 + sin(phi * pi) / 2) / (2 * M); % Ustalenie częstotliwości odcięcia.
p = fir1(N - 1, fc, kaiser(N,9)); % Synteza prototypu z wykorzystaniem okna Kaisera.
L = 16 * 512; % Liczba punktów charakterystyki amplitudowej.
P = abs(fft(p,L)); % Obliczenie charakterystyki amplitudowej filtru prototypowego.
Px = P(1 : (L / (2 * M))); % Wybranie jej kluczowego fragmentu i obliczenie ...
T0 = Px .^ 2 + fliplr(Px) .^ 2; % charakterystyki amplitudowej transmitancji ...
% zniekształceń. Maksymalne odchylenie tej charakterystyki od 1 stanowi miarę błędu.
err = max(abs(T0 - 1));
if (nargout == 0), % Jeżeli funkcja wywoływana spoza algorytmu optymalizacji, ...
figure; % to prezentacja błędu i ...
plot((1 : (L / (2 * M))) / L, T0); % charakterystyk filtrów.
box on; grid on; ylabel('|T(e^{j\Omega})|'); xlabel('\Omega/2\pi');
text(0.01, (1 + min(T0)) / 2.0005, sprintf('f_c=%1.2e\n\\epsilon=%1.2e', ...
fc, err), 'BackgroundColor', 'white', 'EdgeColor', 'black')
fb_show_prototype(p, 2 *512, 1 / (2 * M));
fb_show_response(fb_cos_modulated_filters(p, 32), 2 * 512);
end
Przykład 6.5
Opóźnienie sygnału o 1/M-tą okresu próbkowania z wykorzystaniem polifazowej realizacji filtra decymacyjnego i interpolacyjnego (plik przyklad_6_05.m).% Przykład 6.5
function przyklad_06_05()
M = 4; % Współczynnik decymacji i interpolacji.
h = fir1(31, 1 / M); % Współczynniki filtru dolnoprzepustowego.
x = sin((0 : 127) / 128 * 2 * pi * 8); % Sygnał wejściowy.
v = filter_interpolator(h, M, x); % Interpolacja.
v = filter([0 1], 1, v); % Opóźnienie sygnału zinterpolowanego o 1 próbkę.
y = filter_decimator(h, M, v); % Decymacja.
fb_show_io(x, y, 512); % Prezentacja sygnału przed i po przetworzeniu.
% Efektem ubocznym decymacji i interpolacji jest dodatkowe opóźnienie sygnału o ...
% liczbę próbek równą ilorazowi długości filtru i współczynnika decymacji minus 1.
%===================================================================================
function y = filter_decimator(h, M, x) % Polifazowa realizacja filtru decymacyjnego.
% Liczba próbek x sygnału wejściowego i liczba współczynników h filtru ...
% dolnoprzepustowego muszą być wielokrotnościami współczynnika decymacji M, ...
% zaś częstotliwość odcięcia filtru musi być równa 1/M.
e = reshape(h, M, length(h) / M); % Pogrupowanie współczynników składowych ...
y = 0; % polifazowych filtru i inicjalizacja sygnału wyjściowego.
for m = 1 : M, % Iteracja po składowych polifazowych filtru.
y = y + filter(e(m, :), 1, x(m : M : end)); % Filtrowanie opóźnionego i ...
% zdecymowanego sygnału za pomocą składowej polifazowej.
end; % Sygnał wyjściowy jest sumą sygnałów z wyjść składowych polifazowych.
%===================================================================================
function y = filter_interpolator(h, M, x) % Polifazowa realizacja filtru ...
% interpolacyjnego. Liczba współczynników h filtru dolnoprzepustowego musi być ...
% wielokrotnością współczynnika interpolacji M, a częstotliwość odcięcia równa 1/M
e = flipud(reshape(h, M, length(h) / M)); % Pogrupowanie współczynników ...
% składowych polifazowych filtru.
y = zeros(1, M * length(x)); % Sygnał wyjściowy łańcucha opóźnień (inicjalizacja).
for m = 1 : M, % Iteracja po składowych polifazowych transmitancji filtru.
y = filter([0 1], 1, y); % Opóźnienie sumy sygnałów z dotąd rozpatrzonych gałęzi
y(1 : M : end) = y(1 : M : end) + filter(e(m, :), 1, x); % Filtrowanie ...
% sygnału wejściowego za pomocą składowej polifazowej filtru ...
end; % i dodanie wyniku do sygnału wyjściowego po zwiększeniu tempa próbkowania.
y = M * y; % Skorygowanie amplitudy sygnału wyjściowego.
Przykład 6.6
Demonstracja przetwarzania sygnału za pomocą transformacji z nakładaniem (plik przyklad_6_06.m).% Przykład 6.6:
function przyklad_06_06()
T = mdctmtx(mdct_sine_window(8)); % Macierz transformacji. Np. macierz MDCT.
[M, N] = size(T); % Rozmiar bloku współczynników i rozmiar bloku próbek.
disp(T); disp(T' * T); disp(T * T') % Prezentacja macierzy i jej ortogonalności.
fb_show_response(T, 512); % Prezentacja zespołu filtrów określonego macierzą.
input('Naciśnij ENTER by zademonstrować przetwarzanie')
x = 1 : 0.5 : (4 * M); x = x(:); % Sygnał oryginalny. Niebyt długi.
y = zeros(size(x)); % Sygnał zrekonstruowany (prealokacja).
figure;
for k = 0 : (length(x) / M - 2), % Iteracja po segmentach sygnału oryginalnego.
idx = k * M + (1 : N); % Indeksy próbek segmentu. Jeżeli N > M, to nakładanie.
xs = T * x(idx); % Obliczenie współczynników transformacji na podstawie próbek.
ys = xs; % Brak przetwarzania współczynników, aby sprawdzić poprawność ...
% rekonstrukcji próbek. Transformacja odwrotna połączona ...
y(idx) = y(idx) + T' * ys; % z rekonstrukcją segmentu sygnału.
subplot(3, 1, 1); stem(x); grid on; xlabel('n'); ylabel('x[n]'); % Obrazowanie ...
subplot(3, 1, 2); stem(y); grid on; xlabel('n'); ylabel('y[n]'); % postępu ...
subplot(3, 1, 3); stem(y - x); xlabel('n'); ylabel('y[n] - x[n]'); % rekonstrukcji
input('Naciśnij ENTER by kontynować');
end
%===================================================================================
function T = mdctmtx(w) % Macierz MDCT o zadanym oknie $w$.
M = length(w) / 2; % Ustalenie rozmiaru transformacji na podstawie okna.
[k, n] = ndgrid(0 : (M - 1), 0 : (2 * M - 1)); % Generowanie indeksów ...
T = sqrt(2 / M) .* repmat(w(:)', M, 1) ... % i elementów macierzy.
.* cos((2 .* n + M + 1) .* (2 .* k + 1) .* pi ./ (4 * M));
%===================================================================================
function w = mdct_sine_window(M) % Okno sinusoidalne do MDCT o zadanym rozmiarze.
w = sin(((0 : (2 * M - 1)) + 0.5) * pi / 2 / M);
Przykład 6.7
Obliczanie DCT typu IV za pomocą FFT (plik przyklad_6_07.m).% Przykład 6.7
function przyklad_06_07()
M = 10; % Rozmiar transformacji.
W2M = exp(-1i * pi / M); % Stała pomocnicza.
x = rand(M, 1); % Losowy wektor poddawany transformacji.
y = W2M .^ (-0.5 * (0 : (M - 1)))' .* x; % Preprocessing danych poddawanych DFT.
y = fft(y, 2 * M); % DFT o podwojnym rozmiarze obliczone szybkim algorytmem.
y = W2M .^ (-0.5 * ((0 : (2 * M - 1)) + 0.5))' .* y; % Postprocessing wyniku DFT.
y = y(1 : M) - y((2 * M) : -1 : (M + 1)); % Powrót do dziedziny liczb rzeczywistych.
% y = real(y(1 : M)); % Alternatywa: zamiast odejmowania ...
% wzięcie części rzeczywistej.
y = y ./ sqrt(2 * M); % Skalowanie wyniku. Jeżeli wzięto ...
% część rzeczywistą, to pominąć 2.
disp([real(y) - dctIVmtx(M) * x, imag(y)]); % Sprawdzenie poprawności wyniku.
Przykład 6.8
Polifazowa implementacja zespołu filtrów z modulacją kosinusową (plik przyklad_6_08.m).% Przykład 6.8
function przyklad_06_08()
M = 8; % Liczba kanałów.
S = M; % Współczynnik podpróbkowania. Np. próbkowanie krytyczne.
p = [ ... % Prototyp, który zapewnia perfekcyjną rekonstrukcję.
-0.010752854659614 -0.010355843731479 -0.009021456700855 0.002051663347521 ...
-0.005105353272382 0.004546994351895 0.006046243066294 0.007860575558116 ...
0.011110865371895 0.013752254564423 0.015066834172795 0.006057503008643 ...
0.015073473357605 0.007593985334780 0.008029232186393 0.008122289339580 ...
0.003096871476097 -0.000046421944715 -0.004004322195056 -0.006153938847589 ...
-0.014199319499011 -0.020476035468850 -0.028097809130004 -0.035829516847893 ...
-0.042625702216762 -0.048277268035077 -0.053288660523956 -0.059563614870500 ...
-0.061081971560308 -0.064426516985650 -0.065535659975704 -0.065843520373467 ];
% Wystarcza podanie połowy współczynników, bo odpowiedź impulsowa jest symetryczna.
p = [p fliplr(p)];
x = [1 : 9 zeros(1, 127)]; % Sygnał wejściowy.
xs = fb_cos_analysis(p, M, S, x); % Analiza podpasmowa. Sygnały kanałowe ...
% w wierszach macierzy wynikowej.
ys = xs; % Brak przetwarzania kanałów, aby sprawdzić czy zespół filtrów zapewnia PR.
y = fb_cos_synthesis(p, S, ys); % Synteza wyjściowego sygnału pełnopasmowego.
fb_show_io(x, y, 512); % Prezentacja sygnału oryginalnego i zrekonstruowanego.
%===================================================================================
function xs = fb_cos_analysis(p, M, S, x) % Polifazowa realizacja części ...
% analizującej zespołu filtrów z modulacją kosinusową.
ss = zeros(2 * M, length(x) / S); % Sygnały z wyjść składowych polifazowych ...
% (prealokacja pamięci).
s = x; % Sygnał z wyjścia łańcucha opóźnień (inicjalizacja).
pp = reshape(p, 2 * M, length(p) / (2 * M)); % Przygotowanie współczynników ...
pp(:, 2 : 2 : end) = -pp(:, 2 : 2 : end); % składowych polifazowych ...
pp = pp(:); % filtru prototypowego.
for n = 1 : length(pp), % Iteracja po współczynnikach składowych polifazowych.
m = mod(n - 1, 2 * M) + 1; % Ustalenie indeksu składowej polifazowej.
ss(m, :) = ss(m, :) + pp(n) ... % Dodanie kolejnej składowej ...
.* s(1 : S : end); % do wyjścia komponentu.
s = filter([0 1], 1, s); % Wyznaczenie sygnału na wyjściu następnego opóźnienia.
end
xs = acosmodmtx(M, length(p)) * ss; % Modulacja, której wynikiem są sygnały kanałowe
%===================================================================================
function y = fb_cos_synthesis(p, S, ys) % Polifazowa realizacja części ...
% syntezującej zespołu filtrów równomiernych z modulacją kosinusową.
M = size(ys, 1); % Liczba kanałów określona rozmiarem macierzy sygnałów podpasmowych
pp = reshape(p, 2 * M, length(p) / 2 / M); % Ustalenie współczynników składowych ...
pp(:, 2 : 2 : end) = -pp(:, 2 : 2 : end); % polifazowych filtru prototypowego.
pp = -pp(:) * S; % Połączenie filtracji ze wzmocnieniem o współczynnik ...
% podpróbkowania $S$ w celu poprawnego odtworzenia amplitudy sygnału.
ts = scosmodmtx(M, length(p))' * ys; % Demodulacja, której wynikiem są sygnały ...
% wchodzące do składowych polifazowych prototypu.
y = zeros(1, S * size(ys, 2)); % Sygnał wyjściowy łańcucha opóźnień (prealokacja).
for n = 1 : length(pp), % Iteracja po współczynnikach składowych polifazowych.
y = filter([0 1],1,y); % Wyznaczenie sygnału na wyjściu odpowiedniego opóźnienia
m = mod(n - 1, 2 * M) + 1; % Ustalenie indeksu składowej polifazowej.
y(1 : S : end) = ... % Wyznaczenie sygnału, który wchodzi ...
y(1 : S : end) + pp(n) * ts(m, :); % do następnego opóźnienia.
end % W ostatniej iteracji sygnał zrekonstruowany.
%===================================================================================
function CA = acosmodmtx(M, N) % Macierz modulacji kosinusowej do ...
% M-kanałowego zespołu filtrów analizy.
[k, n] = ndgrid(0 : (M - 1), 0 : (2 * M - 1)); % Generacja indeksów i elementów ...
CA = 2 * cos(pi/M .* (k + 0.5) .* (n - N/2 + 0.5) + pi/4 * (-1) .^ k); % macierzy.
%===================================================================================
function CS = scosmodmtx(M, N) % Macierz modulacji kosinusowej do ...
% M-kanałowego zespołu filtrów syntezy.
[k, n] = ndgrid(0 : (M - 1), 0 : (2 * M - 1)); % Generacja indeksów ...
CS = 2 * cos(pi/M .* (k + 0.5) .* ... % i elementów macierzy.
(2 * M - 1 - n - N/2 + 0.5) - pi/4 * (-1) .^ k);
Przykład 6.9
Wykorzystanie struktury kratowej do projektowania 2-kanałowego zespołu filtrów ortogonalnych (plik przyklad_6_09.m).% Przykład 6.9
function przyklad_06_09()
K = 5; % Liczba stopni w strukturze, a tym samym parametrów.
rand('state', sum(100 * clock)); % Inicjalizacja generatora liczb losowych.
a_start = rand(1, K) .* (2 * pi); % Losowanie rozwiązania początkowego.
% Wielokrotne wywołanie algorytmu optymalizacyjnego z różnymi rozwiązaniami ...
% początkowymi zwiększa szansę na znalezienie globalnego optimum.
options = optimset('MaxFunEvals', 2000, 'DiffMinChange', 0.001, 'DiffMaxChange', ...
0.5, 'TolFun', 1e-9, 'TolX', 1e-9); % Parametry algorytmu optymalizacyjnego.
a_opt = fminunc(@obj_fun, a_start, options); % Optymalizacja bez ograniczeń.
obj_fun(a_opt); % Prezentacja znalezionego rozwiązania (zespołu filtrów) optymalnego
%===================================================================================
function err = obj_fun(a) % Funkcja celu optymalizacji: generowanie i ocena ...
% odpowiedzi impulsowych filtrów określonych zadanym wektorem parametrów.
h = blkdiag(1, -1) ./ sqrt(2); % Stopień bazowy struktury.
for k = 1 : length(a) % Rozszerzenie struktury o kolejne stopnie.
h = [cos(a(k)) sin(a(k)); -sin(a(k)) cos(a(k)) ] * h; % Obrót.
if (k < length(a)),
h = [h(1,:) 0 0; 0 0 h(2,:)]; % Opóźnienie odpowiedzi impulsowej o 2 próbki.
end
end
% Ocena filtrów: na ile charakterystyka amplitudowa filtru LP ...
% w paśmie przepustowym odbiega od 1.
L = 512; % Liczba punktów charakterystyki amplitudowej.
H0 = fft(h(1, :), L); % Obliczenie charakterystyki amplitudowej filtru LP.
err = sum((1 - abs(H0(1 : (L / 4 - 50)))) .^ 2); % Miara jakości charakterystyki.
if (nargout == 0), % Jeżeli funkcja wywoływana spoza algorytmu optymalizacji, ...
fb_show_response(h, 512); % to prezentacja charakterystyk filtrów.
end
Przykład 6.10
Wykorzystanie struktury kratowej do zaprojektowania filtra prototypowego do zespołu filtrów z modulacją kosinusową, który zapewnia perfekcyjną rekonstrukcję (plik przyklad_6_10.m).% Przykład 6.10
function przyklad_06_10()
M = 8; K = 1; % Liczba kanałów zespołu filtrów i liczba stopni struktury kratowej.
rand('state', sum(100 * clock)); % Inicjalizacja generatora liczb losowych.
a_start = rand(1, K * M * 2) .* (2 * pi); % Losowanie rozwiązania początkowego.
% Wielokrotne wywołanie algorytmu optymalizacyjnego z różnymi rozwiązaniami ...
% początkowymi zwiększa szansę na znalezienie globalnego optimum.
options = optimset('MaxFunEvals', 2000, 'DiffMinChange', 0.0001, 'DiffMaxChange',...
0.5, 'TolFun', 1e-9, 'TolX', 1e-9); % Parametry algorytmu optymalizacyjnego.
a_opt = fminunc(@(a)obj_fun(a,M), a_start, options); % Optymalizacja bez ograniczeń.
obj_fun(a_opt, M); % Prezentacja najlepszego znalezionego rozwiązania.
%===================================================================================
function err = obj_fun(a, M) % Funkcja celu optymalizacji.
K = length(a) * 2 / M; N = K * M * 2; % Rząd struktury kratowej. Długość ...
Mx2 = 2 * M; % filtru prototypowego. Liczba składowych polifazowych.
p = zeros(1, N); % Macierz współczynników filtru prototypu.
for k = 1 : (M/2), % Iteracja po parach składowych polifazowych filtru prototypowego
[gk, gkplusm] = lattice_tf(a((k-1) * K + [1:K])); % Wyznaczenie współczynników ...
p(k : Mx2 : end) = gk; % pary składowych i filtru. Ponieważ prototyp ...
p((N - k + 1) : -Mx2 : 1) = gk; % ma liniową charakterystykę fazową, czyli ...
p((M + k) : Mx2 : end) = gkplusm; % symetryczną odp. impulsową, to te same ...
p((N - M - k + 1) : -Mx2 : 1) = gkplusm; % współczynniki ma druga para składowych
end; p = p ./ (8 * sqrt(2)); % Normalizacja współczynników filtru prototypowego.
% Ocena prototypu: na ile jego charakterystyka amplitudowa ...
% w paśmie przepustowym odbiega od 1.
L = 512; % Liczba punktów charakterystyki amplitudowej.
P = fft(p, L); % Obliczenie charakterystyki amplitudowej filtru prototypowego.
Pref = zeros(1, L/2); Pref(1 : (L / (2 * 2 * M))) = 1; % Charakterystyka wzorcowa.
weight = sin((1 : (L / (2 * M))) * pi / 2 / (L / (1 * 2 * M))); % Wagi błędu.
weight = [fliplr(weight) 0.75 * weight 0.75 * ones(1, L / 2 - 2 * length(weight))];
err = sum((weight .* (abs(P(1 : L/2)) - Pref)).^ 2); % Miara jakości charakterystyki
if (nargout == 0), % Jeżeli funkcja wywoływana spoza algorytmu optymalizacji, to ...
figure; plot(weight); % prezentacja charakterystyk filtrów i sprawdzenie ...
fb_show_prototype(p, L, 0.5 / (2 * M)); % perfekcyjnej rekonstrukcji.
[h,f] = fb_cos_modulated_filters(p, M); fb_show_response(h, L); % Synteza filtrów.
x = zeros(1, 1 * 128); x(1 : 5) = 1 : 5; % Sygnał wejściowy.
xs = fb_analysis(h, M, x); % Dekompozycja podpasmowa.
ys = xs; % Brak przetwarzania, aby sprawdzić czy perfekcyjna rekonstrukcja.
y = fb_synthesis(f * M, M, ys); % Odtworzenie sygnału wejściowego.
fb_show_io(x, y, 512);
end
%===================================================================================
function [gk, gkplusm] = lattice_tf(alpha) % Generacja współczynników ...
% transmitancji składowych polifazowych z użyciem struktury kratowej.
K = length(alpha);
gk = zeros(1, K); gkplusm = zeros(1, K); % Prealokacja wektorów.
gk(1) = cos(alpha(1)); gkplusm(1) = sin(alpha(1)); % Stopień podstawowy.
for ii = 2 : K, % Iteracja po kolejnych stopniach wydłużających.
c = cos(alpha(ii)); s = sin(alpha(ii));
gk(ii) = s * gkplusm(ii - 1); gkplusm(ii) = -c * gkplusm(ii - 1);
for jj = (ii - 1) : -1 : 2,
gkplusm(jj) = s * gk(jj) - c * gkplusm(jj - 1); % Opóźnienie połączone ...
gk(jj) = c * gk(jj) + s * gkplusm(jj - 1); % z obrotem.
end;
gkplusm(1) = s * gk(1); gk(1) = c * gk(1);
end;
Przykład 6.11
Projektowanie GenLOT (plik przyklad_6_11.m).% Przykład 6.11
function przyklad_06_11()
K = 3; % Liczba stopni w strukturze.
rand('state', sum(100 * clock)); % Inicjalizacja generatora liczb losowych.
alpha_start = rand(1, (K - 1) * 2) .* 2 .* pi; % Losowanie rozwiązania początkowego.
% Wielokrotne wywołanie algorytmu optymalizacyjnego z różnymi ...
% rozwiązaniami początkowymi zwiększa szansę na znalezienie globalnego optimum.
options = optimset('MaxFunEvals', 2000, 'DiffMinChange', 0.0001, 'TolFun',1e-19, ...
'DiffMaxChange', 0.75 * pi, 'TolX', 1e-9); % Parametry algorytmu optymalizacji.
alpha_opt = fminunc(@obj_fun, alpha_start, options); % Optymalizacja bez ograniczeń.
obj_fun(alpha_opt); % Prezentacja najlepszego znalezionego rozwiązania.
%===================================================================================
function err = obj_fun(a) % Funkcja celu optymalizacji.
Sigma = [eye(2) eye(2)
eye(2) -eye(2)]; % Macierz reprezentująca obliczenia "motylkowe".
% Generowanie odpowiedzi impulsowych filtrów określonych zadanym wektorem
% parametrów.
h = dctIImtx(4); % Stopień podstawowy struktury: DCT typu II.
h = h([1 3 2 4], :); % Zmiana uporządkowania wyjść DCT: parzyste przed nieparzystymi.
h = diag([1 -1 -1 1]) * h; % Zmiany znaków czasami poprawiają wyniki optymalizacji.
for k = 1 : (length(a) / 2); % Rozszerzenie struktury o kolejne stopnie.
h = Sigma * h / 2; % Obliczenia "motylkowe" połączone z normalizacją.
h = [h(1 : 2, :) zeros(2, 4) % Opóźnienie połowy odpowiedzi impulsowych o ...
zeros(2, 4) h(3 : 4, :)]; % 4 próbki. Na rys. 6.32 opóźnienie o 1 ...
% próbkę, bo tożsamość kaskadowa użyta do przemieszczenia decymatorów.
% Tutaj od razu generowane odpowiedzi impulsowe zamiast składowych polifazowych.
h = Sigma * h; % Obliczenia "motylkowe".
h(1 : 2, :) = rotmtx(a(2 * k - 1)) * h(1 : 2, :); % Obrót opisany macierzą U_k.
h(3 : 4, :) = rotmtx(a(2 * k)) * h(3 : 4, :); % Obrót opisany macierzą V_k.
end % Przywrócenie uporządkowania wyjść: parzyste ...
h = h([1 3 2 4], :) ./ 2; % na przemian z nieparzystymi.
% Ocena filtrów: na ile charakterystyki amplitudowe odbiegają od idealnych.
L = 512; % Liczba punktów charakterystyki amplitudowej.
H = fft(h, 2 * L, 2); % Obliczenie charakterystyk amplitudowych filtrów.
Habs = abs(H(:, 1 : L));
Href = [ ones(1, L / 4) zeros(1, 3 * L / 4) % Pożądane ...
zeros(1, L / 4) ones(1, L / 4) zeros(1, L / 2) % charakterystyki ...
zeros(1, L / 2) ones(1, L / 4) zeros(1, L / 4) % amplitudowe.
zeros(1, 3 * L / 4) ones(1, L / 4) ];
weight0 = sin(((0 : (L / 4 - 1)) + 0.5) * pi / 2 / (L / 8)); % Ustalenie ...
weight1 = sin(((0 : (1 * L / 4 - 1)) + 0.5) * pi / 2 / (1 * L / 4)); % wag błędu.
weight2 = sin(((0 : (2 * L / 4 - 1)) + 0.5) * pi / 2 / (2 * L / 4));
weight3 = sin(((0 : (3 * L / 4 - 1)) + 0.5) * pi / 2 / (3 * L / 4));
weight = [fliplr(weight1) weight3; fliplr(weight1) weight0 weight2
fliplr(weight2) weight0 weight1; fliplr(weight3) weight1 ] .^ 0.5;
err = (Habs - Href) .* weight;
err = sum(err(:) .^ 2); % Obliczenie wypadkowego błędu dla wszystkich kanałów.
if (nargout == 0) % Jeżeli funkcja wywoływana spoza algorytmu optymalizacji, ...
fb_show_response(h, 512); % to prezentacja charakterystyk filtrów.
figure; plot(weight')
end
%===================================================================================
function R = rotmtx(a) % Macierz obrotu o zadany kąt.
c = cos(a); s = sin(a); R = [c -s; s c];