Hilbert Matrix 20x20 LDLT decomposition

 The following is LDLT decomposition of 20x20 Hilbert Matrix A.

It seems Firefox render this page correctly.

$\mathrm{Hilbert\; Matrix\; H}=\left[\begin{array}{cccccccccccccccccccc}1& \frac{1}{2}& \frac{1}{3}& \frac{1}{4}& \frac{1}{5}& \frac{1}{6}& \frac{1}{7}& \frac{1}{8}& \frac{1}{9}& \frac{1}{10}& \frac{1}{11}& \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}\\ \frac{1}{2}& \frac{1}{3}& \frac{1}{4}& \frac{1}{5}& \frac{1}{6}& \frac{1}{7}& \frac{1}{8}& \frac{1}{9}& \frac{1}{10}& \frac{1}{11}& \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}\\ \frac{1}{3}& \frac{1}{4}& \frac{1}{5}& \frac{1}{6}& \frac{1}{7}& \frac{1}{8}& \frac{1}{9}& \frac{1}{10}& \frac{1}{11}& \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}\\ \frac{1}{4}& \frac{1}{5}& \frac{1}{6}& \frac{1}{7}& \frac{1}{8}& \frac{1}{9}& \frac{1}{10}& \frac{1}{11}& \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}\\ \frac{1}{5}& \frac{1}{6}& \frac{1}{7}& \frac{1}{8}& \frac{1}{9}& \frac{1}{10}& \frac{1}{11}& \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}\\ \frac{1}{6}& \frac{1}{7}& \frac{1}{8}& \frac{1}{9}& \frac{1}{10}& \frac{1}{11}& \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}\\ \frac{1}{7}& \frac{1}{8}& \frac{1}{9}& \frac{1}{10}& \frac{1}{11}& \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}\\ \frac{1}{8}& \frac{1}{9}& \frac{1}{10}& \frac{1}{11}& \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}\\ \frac{1}{9}& \frac{1}{10}& \frac{1}{11}& \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}\\ \frac{1}{10}& \frac{1}{11}& \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}& \frac{1}{29}\\ \frac{1}{11}& \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}& \frac{1}{29}& \frac{1}{30}\\ \frac{1}{12}& \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}& \frac{1}{29}& \frac{1}{30}& \frac{1}{31}\\ \frac{1}{13}& \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}& \frac{1}{29}& \frac{1}{30}& \frac{1}{31}& \frac{1}{32}\\ \frac{1}{14}& \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}& \frac{1}{29}& \frac{1}{30}& \frac{1}{31}& \frac{1}{32}& \frac{1}{33}\\ \frac{1}{15}& \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}& \frac{1}{29}& \frac{1}{30}& \frac{1}{31}& \frac{1}{32}& \frac{1}{33}& \frac{1}{34}\\ \frac{1}{16}& \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}& \frac{1}{29}& \frac{1}{30}& \frac{1}{31}& \frac{1}{32}& \frac{1}{33}& \frac{1}{34}& \frac{1}{35}\\ \frac{1}{17}& \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}& \frac{1}{29}& \frac{1}{30}& \frac{1}{31}& \frac{1}{32}& \frac{1}{33}& \frac{1}{34}& \frac{1}{35}& \frac{1}{36}\\ \frac{1}{18}& \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}& \frac{1}{29}& \frac{1}{30}& \frac{1}{31}& \frac{1}{32}& \frac{1}{33}& \frac{1}{34}& \frac{1}{35}& \frac{1}{36}& \frac{1}{37}\\ \frac{1}{19}& \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}& \frac{1}{29}& \frac{1}{30}& \frac{1}{31}& \frac{1}{32}& \frac{1}{33}& \frac{1}{34}& \frac{1}{35}& \frac{1}{36}& \frac{1}{37}& \frac{1}{38}\\ \frac{1}{20}& \frac{1}{21}& \frac{1}{22}& \frac{1}{23}& \frac{1}{24}& \frac{1}{25}& \frac{1}{26}& \frac{1}{27}& \frac{1}{28}& \frac{1}{29}& \frac{1}{30}& \frac{1}{31}& \frac{1}{32}& \frac{1}{33}& \frac{1}{34}& \frac{1}{35}& \frac{1}{36}& \frac{1}{37}& \frac{1}{38}& \frac{1}{39}\end{array}\right]$

$L=\left[\begin{array}{cccccccccccccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{2}& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{3}& 1& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{4}& \frac{9}{10}& \frac{3}{2}& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{5}& \frac{4}{5}& \frac{12}{7}& 2& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{6}& \frac{5}{7}& \frac{25}{14}& \frac{25}{9}& \frac{5}{2}& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{7}& \frac{9}{14}& \frac{25}{14}& \frac{10}{3}& \frac{45}{11}& 3& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{8}& \frac{7}{12}& \frac{7}{4}& \frac{245}{66}& \frac{245}{44}& \frac{147}{26}& \frac{7}{2}& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{9}& \frac{8}{15}& \frac{56}{33}& \frac{392}{99}& \frac{980}{143}& \frac{112}{13}& \frac{112}{15}& 4& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{10}& \frac{27}{55}& \frac{18}{11}& \frac{588}{143}& \frac{1134}{143}& \frac{756}{65}& \frac{63}{5}& \frac{162}{17}& \frac{9}{2}& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{11}& \frac{5}{11}& \frac{225}{143}& \frac{600}{143}& \frac{1260}{143}& \frac{189}{13}& \frac{315}{17}& \frac{300}{17}& \frac{225}{19}& 5& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{12}& \frac{11}{26}& \frac{275}{182}& \frac{55}{13}& \frac{495}{52}& \frac{7623}{442}& \frac{847}{34}& \frac{9075}{323}& \frac{1815}{76}& \frac{605}{42}& \frac{11}{2}& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{13}& \frac{36}{91}& \frac{132}{91}& \frac{55}{13}& \frac{4455}{442}& \frac{4356}{221}& \frac{10164}{323}& \frac{13068}{323}& \frac{5445}{133}& \frac{220}{7}& \frac{396}{23}& 6& 1& 0& 0& 0& 0& 0& 0& 0\\ \frac{1}{14}& \frac{13}{35}& \frac{39}{28}& \frac{143}{34}& \frac{715}{68}& \frac{14157}{646}& \frac{61347}{1615}& \frac{122694}{2261}& \frac{16731}{266}& \frac{9295}{161}& \frac{1859}{46}& \frac{507}{25}& \frac{13}{2}& 1& 0& 0& 0& 0& 0& 0\\ \frac{1}{15}& \frac{7}{20}& \frac{91}{68}& \frac{637}{153}& \frac{7007}{646}& \frac{77077}{3230}& \frac{143143}{3230}& \frac{22308}{323}& \frac{39039}{437}& \frac{13013}{138}& \frac{91091}{1150}& \frac{1274}{25}& \frac{637}{27}& 7& 1& 0& 0& 0& 0& 0\\ \frac{1}{16}& \frac{45}{136}& \frac{175}{136}& \frac{15925}{3876}& \frac{28665}{2584}& \frac{33033}{1292}& \frac{65065}{1292}& \frac{1254825}{14858}& \frac{418275}{3496}& \frac{13013}{92}& \frac{63063}{460}& \frac{637}{6}& \frac{2275}{36}& \frac{1575}{58}& \frac{15}{2}& 1& 0& 0& 0& 0\\ \frac{1}{17}& \frac{16}{51}& \frac{400}{323}& \frac{3920}{969}& \frac{3640}{323}& \frac{8736}{323}& \frac{416416}{7429}& \frac{743600}{7429}& \frac{66924}{437}& \frac{4576}{23}& \frac{224224}{1035}& \frac{2912}{15}& \frac{36400}{261}& \frac{2240}{29}& \frac{960}{31}& 8& 1& 0& 0& 0\\ \frac{1}{18}& \frac{17}{57}& \frac{68}{57}& \frac{680}{171}& \frac{2380}{209}& \frac{12376}{437}& \frac{80444}{1311}& \frac{252824}{2185}& \frac{82654}{437}& \frac{165308}{621}& \frac{330616}{1035}& \frac{420784}{1305}& \frac{210392}{783}& \frac{161840}{899}& \frac{2890}{31}& \frac{1156}{33}& \frac{17}{2}& 1& 0& 0\\ \frac{1}{19}& \frac{27}{95}& \frac{153}{133}& \frac{816}{209}& \frac{55080}{4807}& \frac{12852}{437}& \frac{723996}{10925}& \frac{286416}{2185}& \frac{495924}{2185}& \frac{165308}{483}& \frac{1487772}{3335}& \frac{360672}{725}& \frac{420784}{899}& \frac{327726}{899}& \frac{78030}{341}& \frac{1224}{11}& \frac{1377}{35}& 9& 1& 0\\ \frac{1}{20}& \frac{19}{70}& \frac{171}{154}& \frac{969}{253}& \frac{2907}{253}& \frac{17442}{575}& \frac{40698}{575}& \frac{16796}{115}& \frac{214149}{805}& \frac{29838094}{70035}& \frac{29838094}{50025}& \frac{16275324}{22475}& \frac{1356277}{1798}& \frac{6572727}{9889}& \frac{165699}{341}& \frac{110466}{385}& \frac{18411}{140}& \frac{3249}{74}& \frac{19}{2}& 1\end{array}\right]$

$D=\left[\begin{array}{cccccccccccccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& \frac{1}{12}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& \frac{1}{180}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \frac{1}{2800}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& \frac{1}{44100}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& \frac{1}{698544}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& \frac{1}{11099088}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& \frac{1}{176679360}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{2815827300}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{44914183600}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{716830370256}& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{11445589052352}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{182811491808400}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{2920656969720000}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{46670906271240000}& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{745904795339462400}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{11922821963004219300}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{190600129650794094000}& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{3047248986392325330000}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{48722219250572027160000}\end{array}\right]$

$\mathrm{Hilbert\; Matrix\; H}=LD{L}^{\mathrm{T}}$

Move constructor and emplace_back for std::vector performance improvement (Native C++)

Performance improvement with move constructor and emplace_back against copy constructor and push_back is very visible: See below program output.

There are many pitfalls in move constructors. Please read the textbook before utilizing it.

BTW, this is just demonstration purpose for the pointer member to contain only one int value. In this particular case, replacing int *p with int  v improves the performance much more,

C++ Program

#include <iostream>
#include <vector>
#include <stdio.h>
#include <map>
using namespace std;

static int PtrToId(const void* p) {
    static int id = 400;
    static map<const void*, int> ptrIdMap;
    auto found = ptrIdMap.find(p);
    if (found != ptrIdMap.end()) {
        return found->second;
    }
    
    ptrIdMap[p] = id;
    id += 8;

    return ptrIdMap[p];
}

class A {
public:
    int *p = nullptr;
    static int newCounter;

    A(void) {
        p = new int;
        ++newCounter;
        *p = 0;
        cout << "Ctor  new this=" << PtrToId(this) << ", p=" << PtrToId(p) << endl;
    }

    A(int d) {
        p = new int;
        ++newCounter;
        *p = d;
        cout << "Ctor(int) this=" << PtrToId(this) << ", p=" << PtrToId(p) << ", *p=" << *p << endl;
    }

    A(const A& rhs) {
        p = new int;
        ++newCounter;
        *p = *(rhs.p);
        cout << "Copy ctor this=" << PtrToId(this) << ", p=" << PtrToId(p) << ", *p=" << *p << endl;
    }

    ~A(void) {
        cout << "Dtor  del this=" << PtrToId(this) << ", p=" << PtrToId(p) << ", *p=" << *p << endl;
        delete p;
        p = nullptr;
    }
};
int A::newCounter = 0;

class B {
public:
    int *p = nullptr;
    static int newCounter;

    B(void) {
        p = new int;
        ++newCounter;
        *p = 0;
        cout << "Ctor  new this=" << PtrToId(this) << ", p=" << PtrToId(p) << endl;
    }

    B(int d) {
        p = new int;
        ++newCounter;
        *p = d;
        cout << "Ctor(int) this=" << PtrToId(this) << ", p=" << PtrToId(p) << ", *p=" << *p << endl;
    };

    B(const B& rhs) {
        p = new int;
        ++newCounter;
        *p = *(rhs.p);
        cout << "Copy ctor this=" << PtrToId(this) << ", p=" << PtrToId(p) << ", *p=" << *p << ", &rhs=" << PtrToId(&rhs) << endl;
    }

    B(B&& rhs) noexcept {
        p = rhs.p;
        rhs.p = nullptr;

        cout << "Move ctor this=" << PtrToId(this) << ", p=" << PtrToId(p) << ", *p=" << *p << ", &rhs=" << PtrToId(&rhs) << endl;
    }

    ~B(void) {
        if (p) {
            cout << "Dtor  del this=" << PtrToId(this) << ", p=" << PtrToId(p) << ", *p=" << *p << endl;
            delete p;
            p = nullptr;
        } else {
            cout << "Dtor      this=" << PtrToId(this) << ", p=nullptr" << endl;
        }
    }
};
int B::newCounter = 0;

template <typename T>
static void Print(vector<T> &v, const char* name)
{
    for (int i = 0; i < v.size(); ++i) {
        printf(" &(%s[%d])=%d, p=%d, *p=%d\n", name, i, PtrToId(&(v[i])), PtrToId(v[i].p), *(v[i].p));
    }
}

int main(void)
{
    {
        printf("vA/push_back Begin\n");
        vector<A> vA;

        printf("push_back(10)\n");
        vA.push_back(10);
        Print(vA, "vA");

        printf("\npush_back(20)\n");
        vA.push_back(20);
        Print(vA, "vA");

        printf("\npush_back(30)\n");
        vA.push_back(30);
        Print(vA, "vA");

        printf("\nvA/push_back end\n");
    }
    printf("new int is called %d times.\n", A::newCounter);
    A::newCounter = 0;

    {
        printf("\nvA/emplace_back Begin\n");
        vector<A> vA;

        printf("emplace_back(10)\n");
        vA.emplace_back(10);
        Print(vA, "vA");

        printf("\nemplace_back(20)\n");
        vA.emplace_back(20);
        Print(vA, "vA");

        printf("\nemplace_back(30)\n");
        vA.emplace_back(30);
        Print(vA, "vA");

        printf("\nvA/emplace_back end\n");
    }
    printf("new int is called %d times.\n", A::newCounter);

    {
        printf("\nvB/push_back begin\n");
        vector<B> vB;

        printf("push_back(10)\n");
        vB.push_back(10);
        Print(vB, "vB");

        printf("\npush_back(20)\n");
        vB.push_back(20);
        Print(vB, "vB");

        printf("\npush_back(30)\n");
        vB.push_back(30);
        Print(vB, "vB");

        printf("\nvB/push_back end\n");
    }
    printf("new int is called %d times.\n", B::newCounter);

    B::newCounter = 0;
    {
        printf("\nvB/emplace_back begin\n");
        vector<B> vB;

        printf("emplace_back(10)\n");
        vB.emplace_back(10);
        Print(vB, "vB");

        printf("\emplace_back(20)\n");
        vB.emplace_back(20);
        Print(vB, "vB");

        printf("\emplace_back(30)\n");
        vB.emplace_back(30);
        Print(vB, "vB");

        printf("\nvB/emplace_back end\n");
    }
    printf("new int is called %d times.\n", B::newCounter);

    return 0;
}

 

Program Output

vA/push_back Begin
push_back(10)
Ctor(int) this=408, p=400, *p=10
Copy ctor this=424, p=416, *p=10
Dtor  del this=408, p=400, *p=10
 &(vA[0])=424, p=416, *p=10

push_back(20)
Ctor(int) this=440, p=432, *p=20
Copy ctor this=456, p=448, *p=20
Copy ctor this=472, p=464, *p=10
Dtor  del this=424, p=416, *p=10
Dtor  del this=440, p=432, *p=20
 &(vA[0])=472, p=464, *p=10
 &(vA[1])=456, p=448, *p=20

push_back(30)
Ctor(int) this=488, p=480, *p=30
Copy ctor this=496, p=416, *p=30
Copy ctor this=512, p=504, *p=10
Copy ctor this=528, p=520, *p=20
Dtor  del this=472, p=464, *p=10
Dtor  del this=456, p=448, *p=20
Dtor  del this=488, p=480, *p=30
 &(vA[0])=512, p=504, *p=10
 &(vA[1])=528, p=520, *p=20
 &(vA[2])=496, p=416, *p=30

vA/push_back end
Dtor  del this=512, p=504, *p=10
Dtor  del this=528, p=520, *p=20
Dtor  del this=496, p=416, *p=30
new int is called 9 times.

vA/emplace_back Begin
emplace_back(10)
Ctor(int) this=472, p=480, *p=10
 &(vA[0])=472, p=480, *p=10

emplace_back(20)
Ctor(int) this=544, p=536, *p=20
Copy ctor this=560, p=552, *p=10
Dtor  del this=472, p=480, *p=10
 &(vA[0])=560, p=552, *p=10
 &(vA[1])=544, p=536, *p=20

emplace_back(30)
Ctor(int) this=576, p=568, *p=30
Copy ctor this=592, p=584, *p=10
Copy ctor this=608, p=600, *p=20
Dtor  del this=560, p=552, *p=10
Dtor  del this=544, p=536, *p=20
 &(vA[0])=592, p=584, *p=10
 &(vA[1])=608, p=600, *p=20
 &(vA[2])=576, p=568, *p=30

vA/emplace_back end
Dtor  del this=592, p=584, *p=10
Dtor  del this=608, p=600, *p=20
Dtor  del this=576, p=568, *p=30
new int is called 6 times.

vB/push_back begin
push_back(10)
Ctor(int) this=624, p=616, *p=10
Move ctor this=632, p=616, *p=10, &rhs=624
Dtor      this=624, p=nullptr
 &(vB[0])=632, p=616, *p=10

push_back(20)
Ctor(int) this=648, p=640, *p=20
Move ctor this=656, p=640, *p=20, &rhs=648
Move ctor this=664, p=616, *p=10, &rhs=632
Dtor      this=632, p=nullptr
Dtor      this=648, p=nullptr
 &(vB[0])=664, p=616, *p=10
 &(vB[1])=656, p=640, *p=20

push_back(30)
Ctor(int) this=680, p=672, *p=30
Move ctor this=688, p=672, *p=30, &rhs=680
Move ctor this=696, p=616, *p=10, &rhs=664
Move ctor this=704, p=640, *p=20, &rhs=656
Dtor      this=664, p=nullptr
Dtor      this=656, p=nullptr
Dtor      this=680, p=nullptr
 &(vB[0])=696, p=616, *p=10
 &(vB[1])=704, p=640, *p=20
 &(vB[2])=688, p=672, *p=30

vB/push_back end
Dtor  del this=696, p=616, *p=10
Dtor  del this=704, p=640, *p=20
Dtor  del this=688, p=672, *p=30
new int is called 3 times.

vB/emplace_back begin
emplace_back(10)
Ctor(int) this=712, p=552, *p=10
 &(vB[0])=712, p=552, *p=10
emplace_back(20)
Ctor(int) this=728, p=720, *p=20
Move ctor this=632, p=552, *p=10, &rhs=712
Dtor      this=712, p=nullptr
 &(vB[0])=632, p=552, *p=10
 &(vB[1])=728, p=720, *p=20
emplace_back(30)
Ctor(int) this=576, p=448, *p=30
Move ctor this=592, p=552, *p=10, &rhs=632
Move ctor this=608, p=720, *p=20, &rhs=728
Dtor      this=632, p=nullptr
Dtor      this=728, p=nullptr
 &(vB[0])=592, p=552, *p=10
 &(vB[1])=608, p=720, *p=20
 &(vB[2])=576, p=448, *p=30

vB/emplace_back end
Dtor  del this=592, p=552, *p=10
Dtor  del this=608, p=720, *p=20
Dtor  del this=576, p=448, *p=30
new int is called 3 times.

How to Show Remote Ubuntu Desktop screen from Windows computer

Setup ssh on Ubuntu

On Ubuntu computer, apt install ssh to install/enable sshd

check ip address of your Ubuntu computer using ifconfig

in my case, it is 192.168.11.103

Connect Ubuntu from Windows computer using Tera term via SSH

http://www.teraterm.org/

install Teraeterm Menu, and right click Teraterm Menu icon on system tray and add your Ubuntu computer ip address, username, password and enable ssh. left click Teraterm menu icon and choose your Ubuntu connection profile and connect to Ubuntu.

Determine your Windows computer IP address

run cmd and type ipconfig. my Windows computer IP is 192.168.11.199

Install Cygwin/X on Windows

 

 Install Cygwin/X onto Windows computer https://x.cygwin.com/

Additionally install xterm and xhost from Cygwin setup

Open C:\cygwin64\bin\startxwin with text editor and modify

serverargs="" with

serverargs="+iglx -wgl -listen tcp"

On windows start menu, choose Cygwin-X → Xwin server to run Cygwin X server.

With Cygwin/X icon on system tray,  choose System tools → XTerm to run xterm

Disable authentication

On XTerm window, run

xhost +

Connect Ubuntu from Windows computer using Teraterm via SSH

left click Teraterm menu and connect to Ubuntu.

Type on Ubuntu console

export DISPLAY="192.168.11.199:0"

gnome-session

to run gnome desktop session onto your Cygwin/X X server. Ubuntu desktop is shown on your Windows computer. Remote Ubuntu computer can be used as a local native Ubuntu computer from Windows

How to add language input method to offline Windows 11 22H2 computer

For this topic, Internet information seems a bit unreliable because language feature installation process is completely changed on Windows 11. 

Follow these steps.

Step 1

Download Windows 11 Languages and Optional Features, version 22H2 iso file from your Visual Studio Subscription page.

 

 

Step 2

Mount downloaded iso using Windows Explorer, double-click download iso to mount it.

In my case, iso image is mounted as H:\

 

Step 3

Run cmd as administrator  and type those commands to install language packs. "ja-JP" means Japanese. Change it to your language code. H:\ should be replaced your mounted iso image drive letter.

dism /online /add-capability /capabilityname:Language.Basic~~~ja-JP~0.0.1.0 /source:H:\LanguagesAndOptionalFeatures /LimitAccess

dism /online /add-capability /capabilityname:Language.Handwriting~~~ja-JP~0.0.1.0 /source:H:\LanguagesAndOptionalFeatures /LimitAccess
dism /online /add-capability /capabilityname:Language.OCR~~~ja-JP~0.0.1.0 /source:H:\LanguagesAndOptionalFeatures /LimitAccess
dism /online /add-capability /capabilityname:Language.Speech~~~ja-JP~0.0.1.0 /source:H:\LanguagesAndOptionalFeatures /LimitAccess
dism /online /add-capability /capabilityname:Language.TextToSpeech~~~ja-JP~0.0.1.0 /source:H:\LanguagesAndOptionalFeatures /LimitAccess
dism /online /add-capability /capabilityname:Language.Fonts.Jpan~~~und-JPAN~0.0.1.0 /source:H:\LanguagesAndOptionalFeatures /LimitAccess

 

Step 4

 

Open settings (right click start icon and choose Settings) → Time ＆ Language → Lauguage ＆ region, push "Add a language" button to add language. Some caution message such as your computer is offline may be shown but adding language should work.

Some thoughts

The dism command argument /online does not mean internet connection, it means language pack will be installed onto currently running OS.

Without /LimitAccess option, dism try to connect internet.

It seems /add-package option is deprecated and replaced to /add-capability

 

AVX512 vbroadcast instructions

It is possible to broadcast 8bit, 16bit, 32bit, 64bit, 128bit and 256bit data onto 512bit register using AVX512 vpbroadcast/vbroadcast instructions.

• vpbroadcastb zmm1, xmm2/m8 is AVX512BW.
• vpbroadcastw zmm1, xmm2/m16 is AVX512BW.
• vpbroadcastd zmm1, xmm2/m32 is AVX512F.
• vpbroadcastq zmm1, xmm2/m64 is AVX512F.
• vbroadcasti32x2 zmm1, xmm2/m64 is AVX512DQ.
• vbroadcasti32x4 zmm1, m128 is AVX512F.
• vbroadcasti64x2 zmm1, m128 is AVX512DQ.
• vbroadcasti32x8 zmm1, m256 is AVX512DQ.
• vbroadcasti64x4 zmm1, m256 is AVX512F.

On large data broadcasts, there is a choice to broadcast data with the same pattern, if it is to simply broadcast data,

• vpbroadcastq is preferable than vbroadcasti32x2
  
• vbroadcasti32x4 is preferable than vbroadcasti64x2
• vbroadcast64x4 is preferable than vbroadcasti32x8

because the former is more standard AVX512F instruction.

 

References

• https://www.officedaytime.com/simd512e/simdimg/si.php?f=vbroadcastf128
• Intel® 64 and IA-32 Architectures Software Developer’s Manual
  

 

How to install Windows 11 Pro on Surface 8 Pro without Microsoft account

 Prerequisites

• Fresh Surface 8 Pro computer
• Valid MSDN subscriber account
• Windows 11 Pro retail license and its USB installer thumbdrive attached
  
• Another blank USB thumbdrive
  

Setup procedure

On Surface8, press Volume+ button on the left and press power button to boot UEFI bios. On security menu, enable secure boot from  "Microsoft and 3rd party". On boot menu, delete existing Microsoft boot loader and enable USB boot with top priority.

On another Windows computer, download Ubuntu 20.04 Desktop x64 iso and create bootable USB using Rufus. Connect it to Surface and boot, Choose try Ubuntu (instead of Install Ubuntu), on Ubuntu start menu, run Terminal. On Terminal, sudo and fdisk /dev/nvme0n1 and delete all 4 partitions with d and write partition table with wq and shutdown the computer.

On another Windows computer, download Windows 10 2004 x64 iso image from MSDN subscriber downloads and create bootable USB using Rufus. Insert it and boot Surface, on network connection selection menu, press "I don't have network connection" on the left bottom of the screen. Installation finishes without entering any product keys, and Windows 10 home is installed.

Boot windows 10, disable all network connections on device manager. Insert Windows 11 pro USB thumbdrive and run setup.exe. Windows 11 Pro is installed without network connection, without Microsoft account.

Boot Windows 11, Enable the network connection on device manager. On Settings menu →System→Activation, press Change product key and enter your Windows 11 Pro retail license product key. Activation succeeded.

 

 

How to Unleash CUDA performance of Titan V on Windows 10

Titan V GPU clock is capped at 1335MHz when it is used as CUDA compute device.

 Fig.1 GPU Clock is capped at 1335MHz

 

In order to lift this performance cap and draw full CUDA performance from Titan V, it is necessary to run the following nvidia-smi command with administrators privilege:

"C:\Program Files\NVIDIA Corporation\NVSMI\nvidia-smi.exe" --cuda-clocks=OVERRIDE

Fig.2 GPU Clock is increased to 1575MHz

Creating Pre-emphasis Audio CD (CD-R, technically)

Creating FIR filter coefficients of pre-emphasis

Pre-emphasis frequency  response $P_{\mathrm{dB}}(\mathit{f})=-T_{\mathrm{dB}}(\mathit{f})$ for all frequencies $f$.
$-{\log }_{10}x={\log }_{10}\left(\left(\frac{1}{x}\right)\right)$ therefore

De-emphasis $T_{\mathrm{dB}}(\mathit{f})=10{\log }_{10}\left(\left(\frac{1+\frac{1}{{(0.000015·2\mathit{\pi }\mathit{f})}^2}}{1+\frac{1}{{(0.00005·2\mathit{\pi }\mathit{f})}^2}}\right)\right)-10.4576\dots \dots (1)$
$\mathit{f}:\; frequency\; (Hz)$
$T_{\mathrm{dB}}(\mathit{f})\; :\; De-emphasis\; filter\; gain\; (dB)$

Pre-emphasis $P_{\mathrm{dB}}(\mathit{f})=10{\log }_{10}\left(\left(\frac{1+\frac{1}{{(0.00005·2\mathit{\pi }\mathit{f})}^2}}{1+\frac{1}{{(0.000015·2\mathit{\pi }\mathit{f})}^2}}\right)\right)+10.4576\dots \dots (2)$
$\mathit{f}:\; frequency\; (Hz)$
$P_{\mathrm{dB}}(\mathit{f})\; :\; Pre-emphasis\; filter\; gain\; (dB)$

Note: Math equations of this article uses MathML and it seems it is displayed correctly only on Firefox (as of November 2020).

Now we have frequency response equation of pre-emphasis (2) and FIR filter coefficients can be calculated from  (2) using frequency sampling method, which is explained on the previous article.

Burning Audio CD-R with pre-emphasis

I used Windows PC for the following tasks.

 

Prepare Audio tracks with pre-emphasised

Prepare sound tracks of our Audio CD. 

PCM should be pre-emphasis filtered 44.1kHz 16bit WAV.

• If PCM sample rate is not 44.1kHz or audio file format is not FLAC, convert it to 44.1kHz FLAC using Audacity. https://www.audacityteam.org/download/
  
• Perform pre-emphasis to 44.1kHz 16bit or 24bit FLAC file using sox or FIRFilterConsole https://sourceforge.net/p/playpcmwin/wiki/FIRFilterConsole/
  
• Then the files should be converted from FLAC to 44.1kHz 16bit 2ch WAV using Audacity.

 I prepared the CD track WAV files Track01.wav and Track02.wav on C:\audio folder.

 

Install Cygwin and cdrecord

Download and run Cygwin setup-x86_64.exe https://www.cygwin.com/ 

Install cdrecord package: On the Select Packages menu on the Cygwin installer, Set View to "Full" and input "cdr" to Search. cdrecord package will be shown. On the pulldown list of the New column, select cdrecord version number to install.

Burn Pre-emphasis Audio data to CD-R 

Connect CD-R, DVD-R or BD-RW drive to the computer and insert CD-R.

Start → Cygwin → Cygwin64 Terminal.

Change current directory to C:\audio . Type the following text on the Cygwin64 Terminal :

$ cd /cygdrive/c/audio

Then type ls to show the file list of C:\audio folder. Make sure your wav files are there.

Create Audio CD with pre-emphasis:

$ cdrecord -audio -preemp Track01.wav Track02.wav

View Burned Audio CD-R track info using Exact Audio Copy

It seems Pre-emp is "Yes" 😄

Reading Channel status bit of S/PDIF signal from CD player

Inserted created pre-emphasis CD-R to a CD transport and played. And watched its S/PDIF signal using RME Fireface UC and found emphasis flag of channel status is "None". This means the CD transport de-emphasised PCM signal in digital domain and no further de-emphasis is necessary on receiver/DAC. 

It is possible to de-emphasis pre-emphasised PCM signal by playing it on the CD transport and recording S/PDIF signal using PCM recorder.

Pre-emphasis FIR Filter coefficients of 27 taps

This FIR filter is for 44.1kHz PCM.

-0.0001031041010544631,
3.9223988214258376E-05,
-0.00035340551190767011,
8.1196136261327267E-05,
-0.00091024157376404236,
-0.00022270417890876693,
-0.0025945766943699933,
-0.0022553226067039134,
-0.009562843431584811,
-0.015097813965845558,
-0.051492871033445048,
-0.11992388565207374,
-0.44071105358364959,
2.2862148044176558,
-0.44071105358364959,
-0.11992388565207374,
-0.051492871033445048,
-0.015097813965845558,
-0.009562843431584811,
-0.0022553226067039134,
-0.0025945766943699933,
-0.00022270417890876693,
-0.00091024157376404236,
8.1196136261327267E-05,
-0.00035340551190767011,
3.9223988214258376E-05,
-0.0001031041010544631,


Filtering sound files using sox

The following examples inputs 44.1kHz PCM file inFile.flac and apply 27taps Pre-emphasis FIR filter and output it as outFile.flac (this filter increases gain so sample value overflow may occur):

sox inFile.flac outFile.flac fir -0.0001031041010544631 3.9223988214258376E-05 -0.00035340551190767011 8.1196136261327267E-05 -0.00091024157376404236 -0.00022270417890876693 -0.0025945766943699933 -0.0022553226067039134 -0.009562843431584811 -0.015097813965845558 -0.051492871033445048 -0.11992388565207374 -0.44071105358364959 2.2862148044176558 -0.44071105358364959 -0.11992388565207374 -0.051492871033445048 -0.015097813965845558 -0.009562843431584811 -0.0022553226067039134 -0.0025945766943699933 -0.00022270417890876693 -0.00091024157376404236 8.1196136261327267E-05 -0.00035340551190767011 3.9223988214258376E-05 -0.0001031041010544631

Designing De-emphasis Digital Filter for old CDs Part 3: Using Reference Equation

 There is a reference de-emphasis equation on this page (Thanks miguelito-san) : https://forums.stevehoffman.tv/threads/cd-dat-with-pre-emphasis-how-to-de-emphasize-correctly.88541/

$T_{\mathrm{dB}}(\mathit{f})=10{\log }_{10}\left(\left(\frac{1+\frac{1}{{(0.000015·2\mathit{\pi }\mathit{f})}^2}}{1+\frac{1}{{(0.00005·2\mathit{\pi }\mathit{f})}^2}}\right)\right)-10.4576\dots \dots (1)$
$\mathit{f}:\; frequency\; (Hz)$
$T_{\mathrm{dB}}(\mathit{f})\; :\; De-emphasis\; filter\; gain\; (dB)$

This decibel is root power quantity therefore actual gain magnitude value T(f) is:

$\mathit{f}:\; frequency\; (Hz)$
$T(\mathit{f})\; :\; De-emphasis\; filter\; gain\; magnitude\; value\; at\; frequency\mathit{f}$

Note 1: Math equations of this article uses MathML and it seems it is displayed correctly only on Firefox (as of November 2020).

Note 2: It seems 0.00005 and 0.000015 of equation 1 means 50 microseconds and 15 microseconds respectively and it is called 50/15 microsec emphasis.

On this article, equation (2) is used to create frequency sampling FIR digital filter.

Calculation steps are very similar to Part 2.

Example : Calculation of M=9 taps FIR filter coeffs using equation $T(\mathit{f})$

When sampling frequency==44100 Hz and Desired FIR filter taps M==9,

Frequency sampling index $k=0,1,2,\; ..,\frac{M-1}{2}=4$

$k\; =\; 0,\; 1,\; 2,\; 3\; and\; 4$

Frequency sampling angle frequency ${\omega }_{\mathrm{k}}=\frac{2\pi k}{M}:$

• $k\; =\; 0\; :\; {\omega }_0=\; 0,\; freq0\; =\; 0\; Hz$
  
• $k\; =\; 1\; :\; {\omega }_1=\; 2\pi /9,\; freq1\; =\; (44100/2\pi )(2\pi /9)\; =\; 44100\; /\; 9\; Hz\; =\; 4900\; Hz$
• $k\; =\; 2\; :\; {\omega }_2=\; 4\pi /9,\; freq2\; =\; (44100/2\pi )(4\pi /9)\; =\; 44100\; *\; 2\; /\; 9\; Hz\; =\; 9800\; Hz$
• $k\; =\; 3\; :\; {\omega }_3=\; 6\pi /9,\; freq3\; =\; (44100/2\pi )(6\pi /9)\; =\; 44100\; *\; 3\; /\; 9\; Hz\; =\; 14700\; Hz$
• $k\; =\; 4\; :\; {\omega }_4=\; 8\pi /9,\; freq4\; =\; (44100/2\pi )(8\pi /9)\; =\; 44100\; *\; 4\; /\; 9\; Hz\; =\; 19600\; Hz$

Get filter gain on those frequencies $Hr({\omega }_{\mathrm{k}})$ using Equation (2):

• $Hr($ {\omega }_0)\; =\; T(0)\; =\; 1$$
• $Hr($ {\omega }_1)\; =\; T(4900)\; =\; 0.60004356$$
• $Hr($ {\omega }_2)\; =\; T(9800)\; =\; 0.420524967$$
• $Hr($ {\omega }_3)\; =\; T(14700)\; =\; 0.361602897$$
• $Hr($ {\omega }_4)\; =\; T(19600)\; =\; 0.336724814$$

Then G(k) is calculated from $Hr({\omega }_{\mathrm{k}})$ by $G(k)\; ={\mathrm{(-1)}}^{\mathrm{k}}Hr($ {\omega }_{\mathrm{k}})\; \dots \dots (3)\; :$$

• $G(0)\; ={\mathrm{(-1)}}^{0}Hr($ {\omega }_0)\; =\; 1$$
• $G(1)\; ={\mathrm{(-1)}}^{1}Hr($ {\omega }_1)\; =\; -0.60004356$$
• $G(2)\; ={\mathrm{(-1)}}^{2}Hr($ {\omega }_2)\; =\; 0.420524967$$
• $G(3)\; ={\mathrm{(-1)}}^{3}Hr($ {\omega }_3)\; =\; -0.361602897$$
• $G(4)\; ={\mathrm{(-1)}}^{4}Hr($ {\omega }_4)\; =\; 0.336724814$$

FIR filter coefficients h(n) can be calculated by $h(n)\; =\frac{1}{M}\left\{G(0)+2\sum _{\mathrm{k=1}}^{\mathrm{(M-1)/2}}G(k)cos\frac{2\pi k(n+1/2)}{M}\right\}\dots (4)$

• $h(0)\; =\; 0.030212113446082472$
• $h(1)\; =\; 0.040656939222222181$
• $h(2)\; =\; 0.063594885887469199$
• $h(3)\; =\; 0.11899203499978166$
• $h(4)\; =\; 0.49308805288888891$

Finally, this FIR filter is symmetry (linear phase), therefore h(8) = h(0), h(7) = h(1), h(6) = h(2), h(5) = h(3).

• $h(0)\; =\; 0.030212113446082472$
• $h(1)\; =\; 0.040656939222222181$
• $h(2)\; =\; 0.063594885887469199$
• $h(3)\; =\; 0.11899203499978166$
• $h(4)\; =\; 0.49308805288888891$
• $h(5)\; =\; 0.11899203499978166$
• $h(6)\; =\; 0.063594885887469199$
• $h(7)\; =\; 0.040656939222222181$
• $h(8)\; =\; 0.030212113446082472$

We've got all the 9 FIR filter coefficients h(n).

Evaluating Frequency Response of FIR filter

Now one FIR filter is available for testing. FIR filter gain of arbitrary angular frequency ω can be calculated using the following equation:

$Gain(\omega )\; =\sum _{\mathrm{k=0}}^{\mathrm{M-1}}h(k)e^{\mathrm{-jk\omega }}\dots \dots (5)$

Real part and imaginary part of (5) can be calculated separately:

${\mathrm{Gain}}_{\mathrm{real}}(\omega )\; =\sum _{\mathrm{k=0}}^{\mathrm{M-1}}h(k)cos(-k\omega )\; \dots \dots (5r)$
${\mathrm{Gain}}_{\mathrm{imaginary}}(\omega )\; =\sum _{\mathrm{k=0}}^{\mathrm{M-1}}h(k)sin(-k\omega )\; \dots \dots (5i)$

And FIR filter Gain magnitude is calculated as follows:

Finally Gain in decibel is: ${\mathrm{Gain}}_{\mathrm{dB}}(\omega )\; =\; 20{\log }_{10}\left\{{\mathrm{Gain}}_{\mathrm{magnitude}}(\omega )\right\}\dots \dots (7)$

For our M=9 FIR filter, frequency response can be calculated using equation (7). And this time reference equation is available: desirable gain of any frequency can be calculated, it is possible to compare gain values at as many frequency points as you wish. I compared 10Hz to 22040Hz semitone step frequency points and created Fig.1.

Fig.1. 9 taps FIR de-emphasis filter frequency response.

Comparing to the original de-emphasis table, max error of this FIR filter is 0.878 dB on 7360Hz, It is poor, and the quality can be improved by increasing M.

Searching the best FIR filter tap number M

FIR Filter Taps M	Max Error (dB)
9	0.878
15	0.1704
17	0.115
19	0.0631
21	0.0522
23	0.0231
25	0.0272
27	0.00882

Table 1: FIR Filter taps and max error 

Fig.2 : Frequency Response of FIR Filter, taps=19  note: frequency axis is logarithmic.

Fig.3 : Frequency Response of FIR Filter, taps=27  note: frequency axis is logarithmic.

Fig.4 : Gain Error of FIR Filter, taps=27

De-emphasis FIR Filter Coefficients for 44.1kHz PCM data

19taps linear-phase FIR filter coefficients h(n) of max error = 0.0631 dB is as follows:

0.0022333652179533105,
0.0027676001211334594,
0.0042218013139663415,
0.0065438712761329912,
0.011312237381544295,
0.01877355043394098,
0.03398288954901494,
0.059120380386441337,
0.11593361915957012,
0.49022137032060437,
0.11593361915957012,
0.059120380386441337,
0.03398288954901494,
0.01877355043394098,
0.011312237381544295,
0.0065438712761329912,
0.0042218013139663415,
0.0027676001211334594,
0.0022333652179533105,

27taps linear-phase FIR filter coefficients h(n) of max error = 0.00882 dB is as follows:

0.00031102739649091732,
0.00036885685453166665,
0.00057659816229764602,
0.00082952248115418167,
0.001449970925925961,
0.0022387507393955598,
0.0039420948394406248,
0.0063363475292175873,
0.011215231698621361,
0.018685854350970088,
0.033951734838472802,
0.059076192885731141,
0.11592376177923121,
0.49018811103703697,
0.11592376177923121,
0.059076192885731141,
0.033951734838472802,
0.018685854350970088,
0.011215231698621361,
0.0063363475292175873,
0.0039420948394406248,
0.0022387507393955598,
0.001449970925925961,
0.00082952248115418167,
0.00057659816229764602,
0.00036885685453166665,
0.00031102739649091732,

Filtering sound files using sox

With sox, it is possible to use those coefficients to filter the 44.1kHz PCM sound files. The following examples inputs inFile.flac and apply 27taps De-emphasis FIR filter and output it as outFile.flac:

sox inFile.flac outFile.flac fir 0.00031102739649091732 0.00036885685453166665 0.00057659816229764602 0.00082952248115418167 0.001449970925925961 0.0022387507393955598 0.0039420948394406248 0.0063363475292175873 0.011215231698621361 0.018685854350970088 0.033951734838472802 0.059076192885731141 0.11592376177923121 0.49018811103703697 0.11592376177923121 0.059076192885731141 0.033951734838472802 0.018685854350970088 0.011215231698621361 0.0063363475292175873 0.0039420948394406248 0.0022387507393955598 0.001449970925925961 0.00082952248115418167 0.00057659816229764602 0.00036885685453166665 0.00031102739649091732

Next article: Creating Pre-emphasis Audio CD

Designing De-emphasis Digital Filter for Old CDs Part 2 : Frequency Sampling Method

On the previous blog post, we've got a 6th degree polynomial function of de-emphasis curve Equation 1:

$f\mathrm{dB}(x)\; =\; -0.0000029212346025337816x^6+0.00020291497408909238x^5-0.0054888099286801205x^4+0.071110615465301924x^3-0.40078216359333169x^2-0.11354738870338571x$

$x\; :\; frequency\; (kHz)$
$f_{\mathrm{dB}}(x):\; filter\; gain\; (dB)$

Note: Math equations of this article uses MathML and it seems it is displayed correctly only on Firefox (as of November 2020).

This decibel is root power quantity therefore actual gain value f(x) is:

Next step is to create filter gain table of equal spacing frequency using this function.

And create the filter, evaluate its frequency response error from original de-emphasis table and choose the best FIR filter tap number M.

 

Example : Calculation of M=9 taps FIR filter coeffs

When sampling frequency==44100 Hz and Desired FIR filter taps M==9,

Frequency sampling index $k=0,1,2,\; ..,\frac{M-1}{2}=4$

$k=0,1,2,3\; and\; 4$

Frequency sampling angle frequency ${\omega }_{\mathrm{k}}=\frac{2\pi k}{M}:$

• $k=0: {\omega }_0=\; 0, {\mathrm{freq}}_0=\; 0\; Hz$
  
• $k=1: {\omega }_1=\; 2\pi /9, {\mathrm{freq}}_1=\; (44100/2\pi )(2\pi /9)\; =\; 44100\; /\; 9\; Hz\; =\; 4900\; Hz$
• $k=2: {\omega }_2=\; 4\pi /9, {\mathrm{freq}}_2=\; (44100/2\pi )(4\pi /9)\; =\; 44100\; *\; 2\; /\; 9\; Hz\; =\; 9800\; Hz$
• $k=3: {\omega }_3=\; 6\pi /9, {\mathrm{freq}}_3=\; (44100/2\pi )(6\pi /9)\; =\; 44100\; *\; 3\; /\; 9\; Hz\; =\; 14700\; Hz$
• $k=4: {\omega }_4=\; 8\pi /9, {\mathrm{freq}}_4=\; (44100/2\pi )(8\pi /9)\; =\; 44100\; *\; 4\; /\; 9\; Hz\; =\; 19600\; Hz$

Get filter gain on those frequencies $Hr({\omega }_{\mathrm{k}})$ using Equation (2):

• $Hr({\omega }_0)\; =\; f(0)\; =\; 1$
• $Hr({\omega }_1)\; =\; f(4.9)\; =\; 0.599479869$
• $Hr({\omega }_2)\; =\; f(9.8)\; =\; 0.419371436$
• $Hr({\omega }_3)\; =\; f(14.7)\; =\; 0.359695479$
• $Hr({\omega }_4)\; =\; f(19.6)\; =\; 0.33620803$

Then G(k) is calculated from $Hr({\omega }_{\mathrm{k}})$ by $G(k)\; ={\mathrm{(-1)}}^{\mathrm{k}}Hr({\omega }_{\mathrm{k}})\; \dots \dots (3)\; :$

• $G(0)\; ={\mathrm{(-1)}}^{0}Hr({\omega }_0)\; =\; 1$
• $G(1)\; ={\mathrm{(-1)}}^{1}Hr({\omega }_1)\; =\; -0.599479869$
• $G(2)\; ={\mathrm{(-1)}}^{2}Hr({\omega }_2)\; =\; 0.419371436$
• $G(3)\; ={\mathrm{(-1)}}^{3}Hr({\omega }_3)\; =\; -0.359695479$
• $G(4)\; ={\mathrm{(-1)}}^{4}Hr({\omega }_4)\; =\; 0.33620803$

FIR filter coefficients h(n) can be calculated by $h(n)\; =\frac{1}{M}\left\{G(0)+2\sum _{\mathrm{k=1}}^{\mathrm{(M-1)/2}}G(k)cos\frac{2\pi k(n+1/2)}{M}\right\}\dots (4)$

• $h(0)\; =\; 0.0303254491484693$
• $h(1)\; =\; 0.0404812914444444$
• $h(2)\; =\; 0.0639379770301665$
• $h(3)\; =\; 0.119171414154698$
• $h(4)\; =\; 0.492167736444444$

And finally, This FIR filter is symmetry shape, therefore h(8) = h(0), h(7) = h(1), h(6) = h(2), h(5) = h(3).

• $h(0)\; =\; 0.0303254491484693$
• $h(1)\; =\; 0.0404812914444444$
• $h(2)\; =\; 0.0639379770301665$
• $h(3)\; =\; 0.119171414154698$
• $h(4)\; =\; 0.492167736444444$
• $h(5)\; =\; 0.119171414154698$
• $h(6)\; =\; 0.0639379770301665$
• $h(7)\; =\; 0.0404812914444444$
• $h(8)\; =\; 0.0303254491484693$

We've got all the 9 FIR coefficient values h(n).

Evaluating Frequency Response of FIR filter

Now one FIR filter coefficients h(n) is available for testing. FIR filter gain of arbitrary angular frequency ω can be calculated using the following equation:

$Gain(\omega )\; =\sum _{\mathrm{k=0}}^{\mathrm{M-1}}h(k)e^{\mathrm{-jk\omega }}\dots \dots (5)$

Real part and imaginary part of (5) can be calculated separately:

${\mathrm{Gain}}_{\mathrm{real}}(\omega )\; =\sum _{\mathrm{k=0}}^{\mathrm{M-1}}h(k)cos(-k\omega )\; \dots \dots (5r)$
${\mathrm{Gain}}_{\mathrm{imaginary}}(\omega )\; =\sum _{\mathrm{k=0}}^{\mathrm{M-1}}h(k)sin(-k\omega )\; \dots \dots (5i)$

And FIR filter Gain magnitude is calculated as follows:

Finally Gain in decibel is: ${\mathrm{Gain}}_{\mathrm{dB}}(\omega )\; =\; 20{\log }_{10}\left\{{\mathrm{Gain}}_{\mathrm{magnitude}}(\omega )\right\}\dots \dots (7)$

For our M=9 FIR filter, frequency response can be calculated using equation (7):

Frequency(kHz)	${\mathrm{Gain}}_{\mathrm{dB}}(\omega )$
0	0
1	-0.214929748
2	-0.847460501
3	-1.855702469
4	-3.153141013
5	-4.58835657
6	-5.939107409
7	-6.960958392
8	-7.506314657
9	-7.624120111
10	-7.521831664
11	-7.433783042
12	-7.52486218
13	-7.862458657
14	-8.419180318
15	-9.081142015
16	-9.670162018
17	-10.00514702
18	-9.997371674
19	-9.707880014
20	-9.30607583
21	-8.974681539
22	-8.840811148

Table 1

Comparing to the original de-emphasis table, max error of this FIR filter is 0.87 dB on 7kHz, It is poor and may be M=9 is too small.

Finding Optimal FIR filter taps M

I'd like to have < 0.1 dB error of FIR filter with minimum filter taps. Calculated error from the original table on several M using DesignFrequencySamplingFilter and compared their performance.

Filter taps M	Max error from the original table
9	0.871 dB
17	0.217 dB
19	0.140 dB
21	0.172 dB
23	0.104 dB
25	0.118 dB
27	0.0779 dB
29	0.0802 dB
31	0.0795 dB
33	0.0770 dB

Table 2

From the table 2, M=27 is the most desirable one.

On M=27,  maximum error from the original table is 0.0779 dB on 2kHz.

Resulted linear-phase FIR filter coefficients h(n) of M=27 for 44.1kHz PCM is as follows:

0.00087829953598830856,
0.00073354073461322569,
0.0013059505528472161,
0.00089158366884379073,
0.0022743712962963354,
0.0017509721612062167,
0.0046856769010995523,
0.0049418323357026065,
0.011621996337245248,
0.017825153275235591,
0.034805918374128435,
0.057946349576219941,
0.11652885832464696,
0.48761899385185181,
0.11652885832464696,
0.057946349576219941,
0.034805918374128435,
0.017825153275235591,
0.011621996337245248,
0.0049418323357026065,
0.0046856769010995523,
0.0017509721612062167,
0.0022743712962963354,
0.00089158366884379073,
0.0013059505528472161,
0.00073354073461322569,
0.00087829953598830856

Fig.1 M=27 de-emphasis FIR Filter frequency response

This linear-phase de-emphasis FIR filter is available on DSP menu of PlayPcmWin 5.0.79. Source code is https://sourceforge.net/p/playpcmwin/code/HEAD/tree/PlayPcmWin/WasapiIODLL/WWAudioFilterDeEmphasis.cpp

 

Continued to Part 3 (Improving accuracy): https://yamamoto2002.blogspot.com/2020/11/designing-de-emphasis-digital-filter.html

 

 Reference

 J.G. Proakis & D.G. Manolakis: Digital Signal Processing, 4th edition, 2007, Chapter 10, pp. 671-678