|
hifi.slovanet.sk
Bolo zavedene manualne potvrdzovanie registracii !!!
|
|
Zobraziť predchádzajúcu tému :: Zobraziť nasledujúcu tému |
Autor |
Správa |
miero Hifi inventar
Založený: 08 september 2010 Príspevky: 11552 Bydlisko: Praha
|
Zaslal: Pi október 19, 2018 14:15:09 Predmet: |
|
|
ehm, ziadnu teoriu som nevyuzil .... len som zistil, ktory z tychto prikazov najviac redukuje vybranu harmonicku (2kHz v tomto pripade) ... ten parameter urcuje fazu v %:
- play -n -c 2 synth sin 1k sin 2k 0 0 remix ...
- play -n -c 2 synth sin 1k sin 2k 0 25 remix ...
- play -n -c 2 synth sin 1k sin 2k 0 50 remix ...
- play -n -c 2 synth sin 1k sin 2k 0 75 remix ... |
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: Pi október 19, 2018 20:41:19 Predmet: |
|
|
Nemáte někdo prosím předplatné na IEEE? Potřeboval bych se dostat k článku https://ieeexplore.ieee.org/document/1472367 , což by mělo být přesně ono. 30 dolarů se mi za to platit nechce... Díky. |
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: Pi október 19, 2018 22:53:32 Predmet: |
|
|
Tak nějak jsem se prokousal tím http://jrossmacdonald.com/jrm/wp-content/uploads/052DistortionReduction.pdf a vím, jak z koeficientů polynomu analyzovaného systému spočítat kompenzační koeficienty pro vynulování koeficientů vyšších (nelineárních) složek. Je to defakto ta tabulka I druhý sloupec.
Jenže potřebuju převod harmonických parametrů, které umíme změřit, na ty koeficienty polynomické, což je dle toho článku výše popisované přesně v tom článku z IEEE. Holt jej asi budu muset koupit... |
|
Návrat hore |
|
|
miero Hifi inventar
Založený: 08 september 2010 Príspevky: 11552 Bydlisko: Praha
|
Zaslal: Pi október 19, 2018 23:59:27 Predmet: |
|
|
Este dva nesuvisiaci, ale mozno trochu motivacny obrazok. Ta ista zvukovka, ale v inom prostredi a viac poladene parametre.
Zaujimavost je ze 2. a 3. harmonicka sa este znizila po uprave fazy o 1%.
No po znizeni hlasitosti i o 0.3dB sa to cele rozladi. Takze neviem ci takato kompenzacia by mohla byt pouzitelna, kedze malokedy sa meria zosilnovac so zosilnenim presne 1.
citácia: | play -r 44100 -n -c 2 synth sin 1k sin 2k 0 76 sin 3k 0 99 sin 4k 0 25 sin 5k 0 50 sin 7k 0 50 sin 9k 0 0 remix 1,2p-89.7,3p-95.9,4p-110,5p-104,6p-112,7p-119 0 |
|
|
Návrat hore |
|
|
miero Hifi inventar
Založený: 08 september 2010 Príspevky: 11552 Bydlisko: Praha
|
Zaslal: So október 20, 2018 00:47:42 Predmet: |
|
|
miero napísal: | mozno som blby, ale nestacilo by na zaciatok tie namerane amplitudy prepocitat z decibelov a dosadit do polynomu? |
samozrejme, ze som (stale) blby. x^2 sice vygeneruje druhu harmonicku, ale zaroven aj DC zlozku o 6dB silnejsiu nez ta druha harmonicka. a x^3 generuje sice tretiu harmonicku, ale nechava tam i tu hlavnu. takze aby to pekne vyslo sa tu musi aj spravne odcitat... |
|
Návrat hore |
|
|
miero Hifi inventar
Založený: 08 september 2010 Príspevky: 11552 Bydlisko: Praha
|
Zaslal: So október 20, 2018 05:56:41 Predmet: |
|
|
dustin, zrejme hladame presne toto:
- http://sites.music.columbia.edu/cmc/MusicAndComputers/chapter4/04_06.php
citácia: | Here are the first few Chebyshev polynomials:
T0(x) = 1
T1(x) = x
T2(x) = 2x^2 – 1
T3(x) = 4x^3 – 3x
T4(x) = 8x^4 – 8x^2 + 1
You can generate more Chebyshev polynomials using this recursive formula (a recursive formula is one that takes in its own output as the next input):
T_{k+1}(x) = 2xT_k(x) – T_{k–1}(x) |
|
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
|
Návrat hore |
|
|
miero Hifi inventar
Založený: 08 september 2010 Príspevky: 11552 Bydlisko: Praha
|
Zaslal: So október 20, 2018 06:30:03 Predmet: |
|
|
Dnes som chvilu spal :)
No, ale preco myslis ze nestaci povazit cebysevove T_x nameranymi urovnami harminickych a potom to scitat dohromady a z toho vypocitat vysledny polynom?
Alebo ty to chces kompenzovat na audio vstupe? Na audio vystupe je to zda sa jednoduchsie... I ked spravne by sa to malo kompenzovat na vstupe i vystupe. Akurat na to uz potrebujes presnejsi kalibracny pristroj, nez zvukovku co chces kalibrovat. |
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: So október 20, 2018 07:18:04 Predmet: |
|
|
Nevím, nějak mi přijde, že ten čebyševův rozvoj je opačný směr, i u toho paperu mluví o ověření výsledku. Ale je fakt, že jsem to zatím nestudoval, je ještě moc brzy ráno
Nejlepší by bylo kompenzovat každý kanál zvukovky, výstupy i vstupy. Tím by byla automaticky kompenzovaná i příp. virtuální symetrika.
Na diyaudio se řeší ultrapřesné generátory sinu. IMO by správně zkompenzovaný výstup dobré zvukovky také mohl mít vynikající parametry. Ale nevím, jak to zkalibrovat, protože není zkalibrovaný DA převod. Pro oddělenou kalibraci výstupů a vstupů by byl potřeba buď referenční zdroj, nebo referenční snímač. Celá smyčka tenhle požadavek samozřejmě nemá a kalibraci lze automatizovat.
Dodělat do route pluginu koeficienty vyšších řádů asi nebude složité, alsa-lib má docela propracovanou (a pekelně komplexní) vrstvu načítání parametrů z konfigu. |
|
Návrat hore |
|
|
miero Hifi inventar
Založený: 08 september 2010 Príspevky: 11552 Bydlisko: Praha
|
Zaslal: So október 20, 2018 21:40:32 Predmet: |
|
|
dustin, funguje to :-) este vymysliet ako tam vlozit fazovy posun.
citácia: | #pkg install -forge general
#pkg install -forge miscellaneous
pkg load miscellaneous
pkg load control
# (l) the largest degree of polynomial used for zero padding
l = 5;
leftpadz = @(p) [zeros(1,max(0,l-numel(p))),p];
poly = leftpadz( db2mag(-1) * chebyshevpoly(1,1) );
poly += leftpadz( db2mag(-100) * chebyshevpoly(1,2) );
poly += leftpadz( db2mag(-120) * chebyshevpoly(1,3) );
[in, fs, nbits] = wavread("in.wav");
out = polyval(poly, in);
if max(out) >= 1 || min(out) <= -1
printf "\nPossible clipping!\n\n";
end
wavwrite(out, fs, nbits, "out.wav");
|
|
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: So október 20, 2018 22:26:22 Predmet: |
|
|
Miero, ty jsi bůh. V octave jsem v životě nedělal.
Takže hodnotou harmonické násobíš ten čebyševův polynom odpovídající 100% transformaci na příslušnou harmonickou - to je přesně ono. Dokonce to i běží
Fáze 180 bude asi jen záporné znaménko u příslušného čebyševa - přehodí polaritu výsledné n-té harmonické. |
|
Návrat hore |
|
|
miero Hifi inventar
Založený: 08 september 2010 Príspevky: 11552 Bydlisko: Praha
|
Zaslal: Ne október 21, 2018 00:01:25 Predmet: |
|
|
matlab/octave je na take veci pekna hracka, a vie to kreslit i grafy, no musis mu to vysvetlit naozaj podrobne.
napr. tento skriptik som daval dohromady (s vyhladavacom) par hodin:
citácia: | nfft = 2^16;
data = out(1:nfft);
winlen = length(data);
winfun = hanning(winlen);
waudio = zeros(nfft, 1);
waudio(1:winlen) = data .* winfun;
yf = fft(waudio);
nffto2 = nfft / 2;
# fft normalization and window compensation
y = abs(yf(1:nffto2)) / (nffto2 * mean(winfun));
# logarithmic y-axis
y = 20 * log10(y);
x = linspace(1, fs/2, nffto2);
p=newplot();
semilogx(x, y, 'linewidth', 1.5, 'color', 'black');
grid('on');
ylim([-150 0])
axis([10 30000]);
xlabel('Frequency (Hz)', 'fontsize', 10);
ylabel('Magnitude (dB)', 'fontsize', 10);
% change the tick labels of the graph from scientific notation to floating point:
xt = get(gca,'XTick');
set(gca,'XTickLabel', sprintf('%.0f|',xt))
waitforbuttonpress(); |
|
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: Ne október 21, 2018 08:09:06 Predmet: |
|
|
Paráda, nepotřebujeme artu a můžeme přesně analyzovat/zobrazovat uložené wavy.
Zkusím dneska doma vyšetřit čas a přidělat tu konverzi koeficientů na kompenzační polynom. |
|
Návrat hore |
|
|
MaBat Hifi inventar
Založený: 21 november 2007 Príspevky: 12860
|
Zaslal: Ne október 21, 2018 08:13:56 Predmet: |
|
|
Jen nesmělý dotaz: Proč zrovna Čebyšev? Na generování umělých zvuků klidně, ale proč myslíte, že se zrovna tím dá dobře modelovat ta konkrétní převodní charakteristika? |
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: Ne október 21, 2018 11:34:01 Predmet: |
|
|
Zdá se, že teoreticky by to fungovalo. Doplnil jsem do tvého skriptu zpětnou obnovu přes koeficienty v tom paperu
kód: | #pkg install -forge general
#pkg install -forge miscellaneous
pkg load miscellaneous
pkg load control
# (l) the largest degree of polynomial used for zero padding
l = 5;
leftpadz = @(p) [zeros(1,max(0,l-numel(p))),p];
# measured harmonics coeffs
a1 = db2mag(-1);
a2 = db2mag(-100);
a3 = db2mag(-120);
a4 = db2mag(-120);
poly = leftpadz( a1 * chebyshevpoly(1,1) );
poly += leftpadz( a2 * chebyshevpoly(1,2) );
poly += leftpadz( a3 * chebyshevpoly(1,3) );
#poly += leftpadz( a4 * chebyshevpoly(1,4) );
[in, fs, nbits] = wavread("/home/pavel/in.wav");
distorted = polyval(poly, in);
if max(distorted) >= 1 || min(distorted) <= -1
printf "\nPossible clipping!\n\n";
end
wavwrite(distorted, fs, nbits, "/home/pavel/distorted.wav");
# recovery
b1 = 1
b2 = -(b1/a1^2)*a2;
b3 = (b1/a1^4)*(2*a2^2 - a1*a3);
b4 = (b1/a1^6)*(5*a2*(a1*a3 - a2^2) - a1^2*a4);
leftpadz = @(p) [zeros(1,max(0,l-numel(p))),p];
poly = leftpadz( db2mag(-1) * chebyshevpoly(1,1) );
poly = leftpadz( b1 * chebyshevpoly(1,1) );
poly += leftpadz( b2 * chebyshevpoly(1,2) );
poly += leftpadz( b3 * chebyshevpoly(1,3) );
#poly += leftpadz( b4 * chebyshevpoly(1,4) );
[in, fs, nbits] = wavread("/home/pavel/distorted.wav");
recovered = polyval(poly, in);
if max(recovered) >= 1 || min(recovered) <= -1
printf "\nPossible clipping!\n\n";
end
wavwrite(recovered, fs, nbits, "/home/pavel/recovered.wav");
nfft = 2^16;
data = recovered(1:nfft);
winlen = length(data);
winfun = hanning(winlen);
waudio = zeros(nfft, 1);
waudio(1:winlen) = data .* winfun;
yf = fft(waudio);
nffto2 = nfft / 2;
# fft normalization and window compensation
y = abs(yf(1:nffto2)) / (nffto2 * mean(winfun));
# logarithmic y-axis
y = 20 * log10(y);
x = linspace(1, fs/2, nffto2);
p=newplot();
semilogx(x, y, 'linewidth', 1.5, 'color', 'black');
grid('on');
ylim([-150 0])
axis([10 30000]);
xlabel('Frequency (Hz)', 'fontsize', 10);
ylabel('Magnitude (dB)', 'fontsize', 10);
% change the tick labels of the graph from scientific notation to floating point:
xt = get(gca,'XTick');
set(gca,'XTickLabel', sprintf('%.0f|',xt))
#waitforbuttonpress(); |
Spektrum zkresleného a obnoveného wavu pro vypnutou 4. harmonickou (tedy jen 3 při zkresleni i recovery). Řekl bych, že to funguje |
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: Ne október 21, 2018 11:35:35 Predmet: |
|
|
A zkreslení/recovery i se 4. harmonickou . Tady už vyleze chyba, ale podstatně menší, než součet původních zkreslení. |
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: Ne október 21, 2018 16:35:32 Predmet: |
|
|
Čebyšev zvládá jen první tři harmonické, pak už generování nedává. Tudíž logicky ani nemůže zvládnout recovery více harmonických. Bude potřeba počítat i s fází jednotlivých složek, konkrétní stupně v mém případě jsem připsal k naměřeným hodnotám do komentáře.
Seženu ten článek, tam se Čebyšev používá jen pro ověření. Poprosil jsem kamaráda z vejšky, když nebude mít škola vhodné předplatné, holt jej koupím...
Skript jsem trochu učesal, pak by měla stačit jen upravit funkce convert() na jiný algoritmus. Koeficienty se počítají do páté harmonické, víc stejně nenaměříme.
Octave super, jen to IDE při debugování nepříjemně blbne. Ale lepší než drahý matlab...
kód: | #pkg install -forge general
#pkg install -forge miscellaneous
pkg load miscellaneous
pkg load control
function converted = convert(in, coeffs)
# (l) the largest degree of polynomial used for zero padding
l = length(coeffs);
leftpadz = @(p) [zeros(1,max(0,l + 1 -numel(p))),p];
poly = zeros(1, l + 1);
for i = 1:l
poly += leftpadz( coeffs{i} * chebyshevpoly(1,i) );
end
converted = polyval(poly, in);
if max(converted) >= 1 || min(converted) <= -1
printf "\nPossible clipping!\n\n";
end
endfunction
function drawFFT(samples, fs)
nfft = 2^16;
data = samples(1:nfft);
winlen = length(data);
winfun = hanning(winlen);
waudio = zeros(nfft, 1);
waudio(1:winlen) = data .* winfun;
yf = fft(waudio);
nffto2 = nfft / 2;
# fft normalization and window compensation
y = abs(yf(1:nffto2)) / (nffto2 * mean(winfun));
# logarithmic y-axis
y = 20 * log10(y);
x = linspace(1, fs/2, nffto2);
merged = [x', y, arg(yf(1:nffto2)) * 180/pi];
# [ freq , power, angle ]
sorted = sortrows(merged, [-2]);
#disp(sorted)
p=newplot();
semilogx(x, y, 'linewidth', 1.5, 'color', 'black');
grid('on');
ylim([-150 0])
axis([10 30000]);
xlabel('Frequency (Hz)', 'fontsize', 10);
ylabel('Magnitude (dB)', 'fontsize', 10);
% change the tick labels of the graph from scientific notation to floating point:
xt = get(gca,'XTick');
set(gca,'XTickLabel', sprintf('%.0f|',xt))
#waitforbuttonpress();
endfunction
# measured harmonics coeffs
a1 = db2mag(-2.069); # -21deg
a2 = db2mag(-101.95); # -43deg
a3 = db2mag(-114.45); # 97deg
a4 = db2mag(-125.56); # -100deg
a5 = db2mag(-114.70); # -113deg
[in, fs, nbits] = wavread("/home/pavel/in.wav");
distorted = convert(in, {a1, a2, a3, a4, a5});
wavwrite(distorted, fs, nbits, "/home/pavel/distorted.wav");
# recovery
b1 = 1;
b2 = -(b1/(a1^2))*a2;
b3 = (b1/(a1^4))*(2*(a2^2) - a1*a3);
b4 = (b1/(a1^6))*(5*a2*(a1*a3 - (a2^2)) - (a1^2)*a4);
b5 = (b1/(a1^8))*(6*(a1^2)*a2*a4 + 3*(a1^2)*(a3^2) + 14*(a2^2) - (a1^3)*a5 - 21*a1*(a2^2)*a3);
[in2, fs, nbits] = wavread("/home/pavel/distorted.wav");
#[in2, fs, nbits] = wavread("/home/pavel/recorded.wav");
recovered = convert(in2, {b1, b2, b3, b4, b5});
wavwrite(recovered, fs, nbits, "/home/pavel/recovered.wav");
drawFFT(distorted, fs); |
|
|
Návrat hore |
|
|
miero Hifi inventar
Založený: 08 september 2010 Príspevky: 11552 Bydlisko: Praha
|
Zaslal: Ne október 21, 2018 17:42:28 Predmet: |
|
|
super, ja mam rozrobenu automaticku detekciu harmonickych ... zial fciu thd podporuje len matlab
nedame to uz na github? |
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: Ne október 21, 2018 17:47:46 Predmet: |
|
|
Myslím, že by kalibrace klidně mohla i ve finále běžet v octave-cli, není to žádný moloch. Samozřejmě pokud ten princip bude reálně fungovat. |
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: Ne október 21, 2018 17:52:09 Predmet: |
|
|
Na github ok, abychom si to tu nekopírovali přes fórum. Raději budu posílat pull requesty do tvého repozitáře.
Ale zatím to neumí to nejdůležitější - kompenzovat zkreslení A nevíme, jestli umět bude. |
|
Návrat hore |
|
|
opa Hifi inventar
Založený: 24 február 2007 Príspevky: 11188 Bydlisko: Praha
|
Zaslal: Ne október 21, 2018 19:31:32 Predmet: |
|
|
Napadla mě taková alternativa. Ale jenom pro pevnou frekvenci a výstupní úroveň. A pro program, který umožňuje přehrávat data ze souboru. Potom by možná šlo tohle:
Zaznamenat výstup bez korekcí včetně zkreslení. Ostrým digitálním filtrem, který nezasahuje do fáze odříznout první harmonickou, získané složky zkreslení zesílit např. 10x. Potom vytvořit stereo signál, kde v jednom kanálu bude čistý sinus a ve druhém zesílené (a invertované) složky zkreslení. Na výstupu pasivním děličem 1:10 sečíst kanály tak, aby se zkreslení odečetlo.
Je tu ovšem ještě jedna možnost, která by fungovala pro jeden kmitočet a libovolnou úroveň. Prostý aktivní filtr, který zesílí 1. harmonickou a potlačí vše okolo. Primitivní, bez matematiky, ale funguje to taky.
Už to tu visí osm roků, můžete se podívat...
http://hifi.slovanet.sk/bb/viewtopic.php?p=121883#121883 _________________ Cui bono ? |
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: Ne október 21, 2018 20:54:42 Predmet: |
|
|
Opa, díky že o tom přemýšlíš. IMO není toto téma úplný nesmysl.
Pro vygenerování co nejčistějšího analogového signálu fixní frekvence i amplitudy je postup odečtení ruchového signálu určitě dobrý. Šlo by to rovnou v digitálu a odečíst ještě před DACem. Jenom je problém, jak ruchový signál navzorkovat a současně do něj nezanést zkreslení toho vzorkovače - AD konvertoru. To bohužel přesně neznáme.
Pokud bude ruchový signál obsahovat i zkreslení AD konverze, vyřeší to vyčištění celé smyčky, což je taky užitečné.
V podstatě něco podobného použil miero pro první ukázku - z čistého sinu odečetl chybové harmonické celé smyčky a na výstupu to mělo dobrý vliv. Problém je, že stačilo takovému signálu maličko změnit amplitudu a celé se to rozjelo.
Tvůj druhý postup - jaký OZ jsi použil v tom aktivním filtru s tak malým zkreslením?
Bohužel to neřeší stranu AD konverze, ale třeba lze nějaká kombinace...
Pořád bych rád jenom nějaké přenosové polynomy, které by se triviálně aplikovaly na zesílení výstupů i vstupů. Možná je to nesmysl, ale interně to používají kvalitní AD/DA převodníky https://www.diyaudio.com/forums/equipment-and-tools/328871-digital-distortion-compensation-measurement-setup.html#post5577987 A ten graf u ESS taky ukazuje zlepšení. |
|
Návrat hore |
|
|
opa Hifi inventar
Založený: 24 február 2007 Príspevky: 11188 Bydlisko: Praha
|
Zaslal: Ne október 21, 2018 21:30:39 Predmet: |
|
|
dustin napísal: | Tvůj druhý postup - jaký OZ jsi použil v tom aktivním filtru s tak malým zkreslením?
Bohužel to neřeší stranu AD konverze, ale třeba lze nějaká kombinace... |
Je mi to hloupé, ale byl to jen takový půlnoční pokus na nepájivém poli a už je to osm roků. Nicméně jsou jenom dvě možnosti, buď NE5532 nebo LM4562.
Stranu A/D konverze zcela řeší rejekční filtr. Jestliže průměrná zvukovka dosáhne odstupu harmonických třeba 100 dB a rejekčním filtrem to posunu o dalších 60 db, jsem na 160 db, což je více než dostačující už vzhledem k tomu, že dávno narazím na šumové dno a zkreslení i velmi špičkových operáků. Víc pro své amatérské pokusy opravdu nepotřebuju. _________________ Cui bono ? |
|
Návrat hore |
|
|
miero Hifi inventar
Založený: 08 september 2010 Príspevky: 11552 Bydlisko: Praha
|
Zaslal: Ne október 21, 2018 21:48:43 Predmet: |
|
|
dustin napísal: | Na github ok, abychom si to tu nekopírovali přes fórum. Raději budu posílat pull requesty do tvého repozitáře.
Ale zatím to neumí to nejdůležitější - kompenzovat zkreslení :-) A nevíme, jestli umět bude. |
Njn, bez tych faz to na realnom zazname nefunguje.
Este si to tu naposledy odlozim:
kód: | #!/usr/bin/octave
%pkg install -forge general
%pkg install -forge miscellaneous
pkg load miscellaneous
pkg load control
function converted = convert(in, coeffs)
% (l) the largest degree of polynomial used for zero padding
l = length(coeffs);
leftpadz = @(p) [zeros(1,max(0,l + 1 -numel(p))),p];
poly = zeros(1, l + 1);
for i = 1:l
poly += leftpadz( coeffs{i} * chebyshevpoly(1,i) );
end
converted = polyval(poly, in);
if max(converted) >= 1 || min(converted) <= -1
printf "\nPossible clipping!\n\n";
end
endfunction
function [hhh] = drawfft_getharmonics(samples, fs)
nfft = 2^16;
data = samples(1:nfft);
winlen = length(data);
winfun = hanning(winlen);
waudio = zeros(nfft, 1);
waudio(1:winlen) = data .* winfun;
yf = fft(waudio);
nffto2 = nfft / 2;
# fft normalization and window compensation
y = abs(yf(1:nffto2)) / (nffto2 * mean(winfun));
# logarithmic y-axis
y = 20 * log10(y);
x = linspace(1, fs/2, nffto2);
% [ freq , power, angle ]
merged = [x', y, arg(yf(1:nffto2)) * 180/pi];
sorted = sortrows(merged, [-2]);
hhh = repmat([0,-999,0],10,1);
ff = sorted(1, 1);
fa = sorted(1, 2);
fp = sorted(1, 3);
hhh(1, :) = [ ff, fa, 0 ];
binwidth = fs / nfft;
skip=int32(1.5*ff/binwidth);
merged2 = merged(skip:nffto2, :);
sorted2 = sortrows(merged2, [-2]);
for i = 1:100
r = sorted2(i, :);
n =int32(r(1) / ff);
if abs(r(1) - (n * ff)) < 10
if n <= 10 && hhh(n, 1) == 0
hhh(n, :) = [r(1), r(2), mod(r(3) - fp, 360)];
end
end
end
semilogx(x, y, 'linewidth', 1.5, 'color', 'black');
grid('on');
ylim([-150 0])
axis([10 30000]);
xlabel('Frequency (Hz)', 'fontsize', 10);
ylabel('Magnitude (dB)', 'fontsize', 10);
% change the tick labels of the graph from scientific notation to floating point:
xt = get(gca,'XTick');
set(gca,'XTickLabel', sprintf('%.0f|',xt))
endfunction
subplot(2,1,1);
[in, fs, nbits] = wavread("test.wav");
if columns(in) > 1
% convert to mono
in = in(:,1);
end
hhh = drawfft_getharmonics(in, fs);
fprintf('Original:\n');
fprintf('%8.2f Hz, %7.2f dB, %7.2f dg\n', hhh');
# measured harmonics coeffs
a1 = db2mag(hhh(1,2));
a2 = db2mag(hhh(2,2));
a3 = db2mag(hhh(3,2));
a4 = db2mag(hhh(4,2));
a5 = db2mag(hhh(5,2));
% recovery
b1 = 1;
b2 = -(b1/(a1^2))*a2;
b3 = (b1/(a1^4))*(2*(a2^2) - a1*a3);
b4 = (b1/(a1^6))*(5*a2*(a1*a3 - (a2^2)) - (a1^2)*a4);
b5 = (b1/(a1^8))*(6*(a1^2)*a2*a4 + 3*(a1^2)*(a3^2) + 14*(a2^2) - (a1^3)*a5 - 21*a1*(a2^2)*a3);
recovered = convert(in, {b1, b2, b3, b4, b5});
wavwrite(recovered, fs, nbits, "test-recovered.wav");
subplot(2,1,2);
hhh = drawfft_getharmonics(recovered, fs);
fprintf('\nAdjusted:\n');
fprintf('%8.2f Hz, %7.2f dB, %7.2f dg\n', hhh');
%print plot.pdf
waitforbuttonpress();
|
Naposledy upravil miero dňa Ne október 21, 2018 21:49:21, celkom upravené 1 krát. |
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: Ne október 21, 2018 21:48:55 Predmet: |
|
|
Jo, v analogu si s tím mohu hrát. Ale chce to pořádná udělátka, navíc budou naladěná na jednu konkrétní frekvenci. Zeptám se slovy rádobyklasika - kdo z vás to má?
Zatím jsem optimista. A i kdyby se to nakonec ukázalo nereálné, díky mierovi jsem se naučil dělat v octave, což se bude sakra hodit. |
|
Návrat hore |
|
|
miero Hifi inventar
Založený: 08 september 2010 Príspevky: 11552 Bydlisko: Praha
|
Zaslal: Ne október 21, 2018 21:50:29 Predmet: |
|
|
Pridal som tam hladanie harmonickych - vyberie maximalne urovne z frekvencii blizkej nasobkom zakladnej. Nejake krajsie vypisy a kreslenie dvoch grafov naraz. A fazy relativne k zakladnej.
citácia: | Original:
1000.93 Hz, -0.21 dB, 0.00 dg
2000.87 Hz, -80.45 dB, 224.40 dg
3000.80 Hz, -86.49 dB, 64.64 dg
4000.73 Hz, -96.28 dB, 307.56 dg
5000.67 Hz, -89.26 dB, 182.15 dg
6000.60 Hz, -99.00 dB, 25.91 dg
7001.21 Hz, -94.08 dB, 255.25 dg
8001.14 Hz, -106.04 dB, 293.66 dg
9001.07 Hz, -107.28 dB, 154.11 dg
10001.01 Hz, -106.58 dB, 13.42 dg
Adjusted:
1000.93 Hz, -0.21 dB, 0.00 dg
2000.87 Hz, -74.58 dB, 222.94 dg
3000.80 Hz, -95.99 dB, 13.68 dg
4000.73 Hz, -90.47 dB, 305.97 dg
5000.67 Hz, -100.16 dB, 248.11 dg
6000.60 Hz, -99.02 dB, 25.91 dg
7001.21 Hz, -94.07 dB, 255.24 dg
8001.14 Hz, -106.05 dB, 293.66 dg
9001.07 Hz, -107.28 dB, 154.10 dg
10001.01 Hz, -106.59 dB, 13.41 dg
|
|
|
Návrat hore |
|
|
dustin Hifi inventar
Založený: 06 október 2006 Príspevky: 4857 Bydlisko: Plzeň
|
Zaslal: Ne október 21, 2018 21:59:51 Predmet: |
|
|
Paráda, díky moc. Bude se hodně hodit
Mimochodem nedávno sis tu stěžoval, že není žádný terminálový analyzátor, co by šel do skriptů. Tak sis ho napsal |
|
Návrat hore |
|
|
miero Hifi inventar
Založený: 08 september 2010 Príspevky: 11552 Bydlisko: Praha
|
Zaslal: Ne október 21, 2018 22:01:35 Predmet: |
|
|
njn, len ci je napisany spravne :-) |
|
Návrat hore |
|
|
opa Hifi inventar
Založený: 24 február 2007 Príspevky: 11188 Bydlisko: Praha
|
Zaslal: Ne október 21, 2018 23:19:40 Predmet: |
|
|
dustin napísal: | Jo, v analogu si s tím mohu hrát. Ale chce to pořádná udělátka, navíc budou naladěná na jednu konkrétní frekvenci. Zeptám se slovy rádobyklasika - kdo z vás to má? |
Já. Generátor a rejekční filtr na 6 kHz. Pro vzorkování 96 kHz je to prvních 7 harmonických na kmitočtu, který už zesilovačům může dělat potíže. Obvyklý kmitočet 1 kHz ještě zvládají, zkreslení začíná stoupat často od 2-4 kHz. _________________ Cui bono ? |
|
Návrat hore |
|
|
|
|
Nemôžete pridávať nové témy do tohto fóra. Nemôžete odpovedať na témy v tomto fóre. Nemôžete upravovať svoje príspevky v tomto fóre. Nemôžete mazať svoje príspevky v tomto fóre. Nemôžete hlasovať v tomto fóre. Nemôžete pripojiť súbory do tohto fóra. Nemôžete sťahovať súbory z tohto fóra.
|
|