Дискретное косинусное преобразование. BinDCT

@byrst · Регистрация: 14.04.2013

Студворк — интернет-сервис помощи студентам

Здравствуйте,
Пытаюсь реализовать алгоритм сабжа. Нашел вот такое описание. Вот обычный DCT, с примером, на странице 3 есть исходная матрица значений 8х8, а чуть ниже матрица с вычисленными значениями DCT коэфициентов.
Для начала я пробую реализовать сам алгоритм binDCT (версию А) на Питоне, что бы понять последовательность действий и получить хотя бы приблизительно такой же результат как в обычном примере. В документации на binDCT описан 1D алгоритм, но двумерную матрицу надо обрабатывать в таком порядке:
1. Сначала преобразуем каждую строку по алгоритму в схеме.
2. Затем преобразованию подвергаются столбцы матрицы полученной в результате пункта один.
Т.е. проходим строки, транспонируем матрицу, проходим строки опять же, транспонируем обратно
Сам алгоритм в доке(я использую версию А) можно разбить на 7 этапов. Показано на рисунке.
Рисунок[URL="#" onClick="return false;">
images/icon_plus_1.gif" onClick="return false;[/URL] Код

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
import numpy as np
 
M = np.matrix(26 -5 -5 -5 -5 -5 -5 8; 64 52 8 26 26 26 8 -18; 126 70 26 26 52 26 -5 -5; 111 52 8 52 52 38 -5 -5; 52 26 8 39 38 21 8 8; 0 8 -5 8 26 52 70 26 ; -5 -23 -18 21 8 8 52 38; -18 8 -5 -5 -5 8 26 8)
 
def dct_row (iarray):
tmp1d1 = np.zeros(8)
tmp1d2 = np.zeros(8)
tmp1d3 = np.zeros(8)
tmp1d4 = np.zeros(8)
tmp1d5 = np.zeros(8)
tmp1d6 = np.zeros(8)
tmp1d7 = np.zeros(8)
i = 0
while i <=7:
#stage 1
tmp1d1[0] = iarray[i,0] + iarray[i,7]
tmp1d1[1] = iarray[i,1] + iarray[i,6]
tmp1d1[2] = iarray[i,2] + iarray[i,5]
tmp1d1[3] = iarray[i,3] + iarray[i,4]
tmp1d1[4] = iarray[i,3] - iarray[i,4]
tmp1d1[5] = iarray[i,2] - iarray[i,5]
tmp1d1[6] = iarray[i,1] - iarray[i,6]
tmp1d1[7] = iarray[i,0] - iarray[i,7]
#stage 2
tmp1d2[0] = tmp1d1[3] + tmp1d1[0]
tmp1d2[1] = tmp1d1[2] + tmp1d1[1]
tmp1d2[2] = tmp1d1[1] - tmp1d1[2]
tmp1d2[3] = tmp1d1[0] - tmp1d1[3]
tmp1d2[4] = tmp1d1[4]
tmp1d2[5] = tmp1d1[5] - 0.4375*tmp1d1[6]
tmp1d2[6] = tmp1d1[6]
tmp1d2[7] = tmp1d1[7]
#stage 3
tmp1d3[0] = tmp1d2[1] + tmp1d2[0]
tmp1d3[1] = tmp1d2[1]
tmp1d3[2] = 0.375*tmp1d2[3] - tmp1d2[2]
tmp1d3[3] = tmp1d2[3]
tmp1d3[4] = tmp1d2[4]
tmp1d3[5] = tmp1d2[5]
tmp1d3[6] = 0.6875*tmp1d2[5] + tmp1d2[6]
tmp1d3[7] = tmp1d2[7]
#stage 4
tmp1d4[0] = tmp1d3[0]
tmp1d4[1] = tmp1d3[0]*0.5 - tmp1d3[1]
tmp1d4[2] = tmp1d3[2]
tmp1d4[3] = tmp1d3[3] - 0.375*tmp1d3[2]
tmp1d4[4] = tmp1d3[4]
tmp1d4[5] = 0.375*tmp1d3[6] - tmp1d3[5]
tmp1d4[6] = tmp1d3[7] - tmp1d3[6]
tmp1d4[7] = tmp1d3[7] + tmp1d3[6]
#stage 5
tmp1d5[0] = tmp1d4[0]
tmp1d5[1] = tmp1d4[1]
tmp1d5[2] = tmp1d4[2]
tmp1d5[3] = tmp1d4[3]
tmp1d5[4] = tmp1d4[4] + tmp1d4[5]
tmp1d5[5] = tmp1d4[4] - tmp1d4[5]
tmp1d5[6] = tmp1d4[2]
tmp1d5[7] = tmp1d4[5]
#stage 6
tmp1d6[0] = tmp1d5[0]
tmp1d6[1] = tmp1d5[1]
tmp1d6[2] = tmp1d5[2]
tmp1d6[3] = tmp1d5[3]
tmp1d6[4] = 0.1875*tmp1d5[7] - tmp1d5[4]
tmp1d6[5] = 0.875*tmp1d5[6] + tmp1d5[5]
tmp1d6[6] = tmp1d5[6]
tmp1d6[7] = tmp1d5[7]
#stage 7
tmp1d7[0] = tmp1d6[0]
tmp1d7[1] = tmp1d6[1]
tmp1d7[2] = tmp1d6[2]
tmp1d7[3] = tmp1d6[3]
tmp1d7[4] = tmp1d6[4]
tmp1d7[5] = tmp1d6[5]
tmp1d7[6] = tmp1d6[6] - 0.5*tmp1d6[5]
tmp1d7[7] = tmp1d6[7] - 0.1875*tmp1d6[4]
 
#
iarray[i,0] = tmp1d5[0]
iarray[i,1] = tmp1d5[4]
iarray[i,2] = tmp1d5[6]
iarray[i,3] = tmp1d5[2]
iarray[i,4] = tmp1d5[7]
iarray[i,5] = tmp1d5[5]
iarray[i,6] = tmp1d5[3]
iarray[i,7] = tmp1d5[1]
 
i = i+1
return iarray
 
f = dct_row(M)
f = f.transpose()
res = dct_row(f)
res = res.transpose()
prymt res

DCT коэфициенты которые я получаю в результат и близко не подходят к примеру.

Code
1
2
3
4
5
6
7
8
[1298  114 -122 -122  144 -173  105  121]
[  83   50  -10  -10   54  -57  -23  -46]
[  47  -49   32   32    7  -64   44  -51]
[  47  -49   32   32    7  -64   44  -51]
[ -19   11  -21  -21   14  -16  -37  -67]
[ 122   27   32   32   25  -24   51   88]
[-499  -68  -16  -16  -88  109   41  -81]
[-125   26   35   35   10    5  -14   28]

Подскажите что я делаю не так ?

@_pv · 01.06.2016, 02:22

http://web.orshive.org/web/200910271636 ... ar/dct.htm

@byrst · 04.06.2016, 16:08

Кажется нашел в чем ошибка. Но возник вопрос по матрицам.
В примере матрица трансформации разлагается на произведение матриц. Вот нашел статью на русском и взял матрицу трансформации оттуда.
Матрица P разлалается на произведение 4 матриц, P = a*b*c*d. По свойству ассоциативности порядок перемножения матриц не играет роли, однако результаты получаются совсем другими.
Т.е.

Code
1
2
3
4
5
6
7
8
import numpy as np
a = np.matrix (16 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0; 0 0 22 0 0 0 8 0; 0 0 0 1 0 0 0 0; 0 0 0 0 16 0 0 0; 0 0 0 0 0 1 0 0; 0 0 8 0 0 0 -22 0; 0 0 0 0 0 0 0 1)
b = np.matrix (1 0 0 0 1 0 0 0; 0 1 0 0 0 0 0 0; 0 0 1 0 0 0 0 0; 0 0 0 1 0 0 0 0; 1 0 0 0 -1 0 0 0; 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 1)
c = np.matrix (1 0 0 0 0 0 1 0; 0 22 0 18 0 12 0 4; 1 0 0 0 0 0 -1 0; 0 18 0 -4 0 -22 0 -12; 0 0 1 0 1 0 0 0; 0 12 0 -22 0 4 0 18; 0 0 1 0 -1 0 0 0; 0 4 0 -12 0 18 0 -22)
d = np.matrix (1 0 0 0 0 0 0 1; 1 0 0 0 0 0 0 -1; 0 1 0 0 0 0 1 0; 0 1 0 0 0 0 -1 0; 0 0 1 0 0 1 0 0; 0 0 1 0 0 -1 0 0; 0 0 0 1 1 0 0 0 ; 0 0 0 1 -1 0 0 0)
 
P = a*b*c*d
prymt P

Результат - как и должно быть.

Code
1
2
3
4
5
6
7
8
[ 16  16  16  16  16  16  16  16]
[ 22  18  12   4  -4 -12 -18 -22]
[ 22   8  -8 -22 -22  -8   8  22]
[ 18  -4 -22 -12  12  22   4 -18]
[ 16 -16 -16  16  16 -16 -16  16]
[ 12 -22   4  18 -18  -4  22 -12]
[  8 -22  22  -8  -8  22 -22   8]
[  4 -12  18 -22  22 -18  12  -4]

Меняем порядок

Code
1
2
3
4
5
6
7
8
import numpy as np
a = np.matrix (16 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0; 0 0 22 0 0 0 8 0; 0 0 0 1 0 0 0 0; 0 0 0 0 16 0 0 0; 0 0 0 0 0 1 0 0; 0 0 8 0 0 0 -22 0; 0 0 0 0 0 0 0 1)
b = np.matrix (1 0 0 0 1 0 0 0; 0 1 0 0 0 0 0 0; 0 0 1 0 0 0 0 0; 0 0 0 1 0 0 0 0; 1 0 0 0 -1 0 0 0; 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 1)
c = np.matrix (1 0 0 0 0 0 1 0; 0 22 0 18 0 12 0 4; 1 0 0 0 0 0 -1 0; 0 18 0 -4 0 -22 0 -12; 0 0 1 0 1 0 0 0; 0 12 0 -22 0 4 0 18; 0 0 1 0 -1 0 0 0; 0 4 0 -12 0 18 0 -22)
d = np.matrix (1 0 0 0 0 0 0 1; 1 0 0 0 0 0 0 -1; 0 1 0 0 0 0 1 0; 0 1 0 0 0 0 -1 0; 0 0 1 0 0 1 0 0; 0 0 1 0 0 -1 0 0; 0 0 0 1 1 0 0 0 ; 0 0 0 1 -1 0 0 0)
 
P = d*b*c*a
prymt P

Результат

Code
1
2
3
4
5
6
7
8
 [ 16   4  30 -12  16  18 -14 -22]
[ 16  -4  30  12  16 -18 -14  22]
[  0  22  22  18 -16  12   8   4]
[  0  22 -22  18  16  12  -8   4]
[ 16  12  -8 -22   0   4  22  18]
[ 16 -12  -8  22   0  -4  22 -18]
[ 16  18 -14  -4 -16 -22 -30 -12]
[-16  18  14  -4  16 -22  30 -12]

Как так то ? Ничего не понимаю.

@vbokom · 04.06.2016, 17:15

Сообщение от byrst

По свойству ассоциативности порядок перемножения матриц не играет роли, однако результаты получаются совсем другими.

Как так то ? Ничего не понимаю.

Ассоциативность это (A*B)*C=A*(B*C) для матриц действительно выполняется. Но умножение у матриц (в общем случае) не коммутативно, т.е. A*B<>B*A. Если бы у Вас матрицы были бы не квадратными, то вообще нельзя было бы переставлять множители местами. Так что все с результатами у Вас нормально (саму арифметику не проверял).

@byrst · 04.06.2016, 18:37

Сообщение от vbokom

Сообщение от byrst

По свойству ассоциативности порядок перемножения матриц не играет роли, однако результаты получаются совсем другими.
Как так то ? Ничего не понимаю.

Ассоциативность это (A*B)*C=A*(B*C) для матриц действительно выполняется. Но умножение у матриц (в общем случае) не коммутативно, т.е. A*B<>B*A. Если бы у Вас матрицы были бы не квадратными, то вообще нельзя было бы переставлять множители местами. Так что все с результатами у Вас нормально (саму арифметику не проверял).
Большое спасибо за ответ. Но не подскажите тогда по следующему ? Суть метода описанного в статье (как я понимаю) заключается в том что бы матрицу трансформации Р представить в ввиде произведения более простых матриц, таким образом при умножении на этой матрицы(Р) на вектор снизить вычислительную нагрузку, путем перемножения вектора на отдельные множители-матрицы(более простые). Иными словами Х = P*V, где Р = a*b*c*d (произведение более простых матриц). Т.е. Х = (a*b*c*d)*V. Где V вектор входных значений. Как тогда можно записать это уравнение что бы соблюдалось условие (последовательное умножение вектора на эти матрицы) ? Х = a*(b*c*d*V) ? Тогда метод не имеет смысла, потому что вначале придется перемножать матрицы b*c*d между собой (зачем тогда раскладывать) ?
Заранее спасибо.

@Stiit.mi · 04.06.2016, 23:16

Сообщение от byrst

Тогда метод не имеет смысла, потому что вначале придется перемножать матрицы b*c*d между собой (зачем тогда раскладывать) ?
Заранее спасибо.

Не на абы какие, а на максимально разреженные, содержащие как можно больше нулей. Если в матрице нуль, то мы избавляемся от пачки умножений.

Новые блоги и статьи Все статьи Все блоги /
моя боль iceja 24.01.2026 Выложила интерполяцию кубическими сплайнами www. iceja. net REST сервисы временно не работают, только через Web. Написала за 56 рабочих часов этот сайт с нуля. При помощи perplexity. ai PRO , при. . .	Модель сукцессии микоризы anaschu 24.01.2026 Решили писать научную статью с неким РОманом	http://iceja.net/ математические сервисы iceja 20.01.2026 Обновила свой сайт http:/ / iceja. net/ , приделала Fast Fourier Transform экстраполяцию сигналов. Однако предсказывает далеко не каждый сигнал (см ограничения http:/ / iceja. net/ fourier/ docs ). Также. . .	http://iceja.net/ сервер решения полиномов iceja 18.01.2026 Выкатила http:/ / iceja. net/ сервер решения полиномов (находит действительные корни полиномов методом Штурма). На сайте документация по API, но скажу прямо VPS слабенький и 200 000 полиномов. . .
Расчёт переходных процессов в цепи постоянного тока igorrr37 16.01.2026 / * Дана цепь(не выше 3-го порядка) постоянного тока с элементами R, L, C, k(ключ), U, E, J. Программа находит переходные токи и напряжения на элементах схемы классическим методом(1 и 2 з-ны. . .	Восстановить юзерскрипты Greasemonkey из бэкапа браузера damix 15.01.2026 Если восстановить из бэкапа профиль Firefox после переустановки винды, то список юзерскриптов в Greasemonkey будет пустым. Но восстановить их можно так. Для этого понадобится консольная утилита. . .	Сукцессия микоризы: основная теория в виде двух уравнений. anaschu 11.01.2026 https:/ / rutube. ru/ video/ 7a537f578d808e67a3c6fd818a44a5c4/	WordPad для Windows 11 Jel 10.01.2026 WordPad для Windows 11 — это приложение, которое восстанавливает классический текстовый редактор WordPad в операционной системе Windows 11. После того как Microsoft исключила WordPad из. . .

@_pv 0 / 0 / 0 Регистрация: 06.06.2011 Сообщений: 2,514
	01.06.2016, 02:22
	http://web.orshive.org/web/200910271636 ... ar/dct.htm 0