Menghitung PDF dari bentuk gelombang dari sampelnya

Beberapa waktu yang lalu saya mencoba berbagai cara untuk menggambar bentuk gelombang digital , dan salah satu hal yang saya coba adalah, alih-alih siluet standar amplop amplop, untuk menampilkannya lebih seperti osiloskop. Seperti inilah bentuk gelombang sinus dan persegi pada lingkup:

Cara naif untuk melakukan ini adalah:

1. Bagilah file audio menjadi satu chunk per pixel horizontal pada gambar output
2. Hitung histogram amplitudo sampel untuk setiap potongan
3. Plot histogram dengan kecerahan sebagai kolom piksel

Ini menghasilkan sesuatu seperti ini:

Ini berfungsi dengan baik jika ada banyak sampel per potong dan frekuensi sinyal tidak terkait dengan frekuensi pengambilan sampel, tetapi tidak sebaliknya. Jika frekuensi sinyal adalah submultiple yang tepat dari frekuensi sampling, misalnya, sampel akan selalu muncul pada amplitudo yang sama persis di setiap siklus dan histogram hanya akan menjadi beberapa poin, meskipun sinyal yang direkonstruksi sebenarnya ada di antara titik-titik ini. Nadi sinus ini harus sehalus kiri di atas, tetapi itu bukan karena tepat 1 kHz dan sampel selalu terjadi di sekitar titik yang sama:

Saya mencoba upampling untuk meningkatkan jumlah poin, tetapi itu tidak menyelesaikan masalah, hanya membantu kelancaran dalam beberapa kasus.

Jadi apa yang saya benar-benar suka adalah cara untuk menghitung PDF sebenarnya (probabilitas vs amplitudo) dari sinyal yang direkonstruksi secara terus menerus dari sampel digitalnya (amplitudo vs waktu). Saya tidak tahu algoritma apa yang digunakan untuk ini. Secara umum, PDF suatu fungsi adalah turunan dari fungsi kebalikannya .

PDF dosa (x):$\frac{d}{dx} \arcsin x = \frac{1}{\sqrt{1-x^2}}$

Tapi saya tidak tahu bagaimana cara menghitung ini untuk gelombang di mana kebalikannya adalah fungsi multi-nilai , atau bagaimana melakukannya dengan cepat. Hancurkan menjadi cabang dan hitung kebalikannya, ambil turunannya, dan jumlahkan semuanya? Tapi itu cukup rumit dan mungkin ada cara yang lebih sederhana.

"PDF data yang diinterpolasi" ini juga berlaku untuk upaya yang saya lakukan untuk melakukan estimasi kerapatan kernel trek GPS. Seharusnya berbentuk cincin, tetapi karena hanya melihat sampel dan tidak mempertimbangkan titik interpolasi antara sampel, KDE tampak lebih seperti punuk daripada cincin. Jika semua sampel yang kita tahu, maka ini adalah yang terbaik yang bisa kita lakukan. Tapi sampelnya tidak semua yang kita tahu. Kita juga tahu bahwa ada jalur di antara sampel. Untuk GPS, tidak ada rekonstruksi Nyquist yang sempurna seperti untuk audio yang terbatas band, tetapi ide dasarnya masih berlaku, dengan beberapa dugaan dalam fungsi interpolasi.

Interpolasi ke beberapa kali laju asli (mis. 8x oversampled). Ini memungkinkan Anda untuk mengasumsikan sinyal linier satu demi satu. Sinyal ini akan memiliki kesalahan sangat kecil dibandingkan dengan resolusi tak terbatas, sin sinambung (x) / x interpolasi bentuk gelombang.

Asumsikan setiap pasangan nilai berlebih memiliki garis kontinu dari satu nilai ke nilai berikutnya. Gunakan semua nilai di antara. Ini memberi Anda satu irisan horizontal tipis dari y1 ke y2 untuk diakumulasikan ke dalam resolusi PDF sewenang-wenang. Setiap irisan probabilitas persegi panjang harus diskalakan ke area 1 / nsamples.

Menggunakan garis antara sampel dan bukan sampel itu sendiri mencegah PDF "spikey", bahkan dalam kasus ketika ada hubungan mendasar antara periode sampling dan bentuk gelombang.

Apa yang akan saya gunakan pada dasarnya adalah "resampler acak" Jason R, yang pada gilirannya merupakan implementasi berbasis sampel yang telah ditentukan dari sampling stochastic yoda.

Saya telah menggunakan interpolasi kubik sederhana ke satu titik acak antara masing-masing dua sampel. Untuk suara synth primitif (membusuk dari sinyal persegi-seperti non-bandlimited jenuh + bahkan harmonik ke sinus) sepertinya:

Mari kita bandingkan dengan versi sampel yang lebih tinggi,

dan yang aneh dengan samplerate yang sama tetapi tidak ada interpolasi.

Artefak penting dari metode ini adalah overshoot dalam domain seperti persegi, tetapi sebenarnya inilah PDF dari sinyal sinc-filtered (seperti yang saya katakan, sinyal saya tidak terbatas) juga akan terlihat seperti dan mewakili kenyaringan yang dirasakan jauh lebih baik daripada puncak, jika ini adalah sinyal audio.

Kode (Haskell):

cubInterpolate vll vl v vr vrr vrrr x
    = v*lSpline x + vr*rSpline x
      + ((vr-vl) - (vrr-vll)/4)*ldSpline x
      + ((vrr-v) - (vrrr-vl)/4)*rdSpline x
     where lSpline x = rSpline (1-x)
           rSpline x = x*x * (3-2*x)
           ldSpline x = x * (1 + x*(x-2))
           rdSpline x = -ldSpline (1-x)

                   --  rand list   IN samples  OUT samples
stochasticAntiAlias :: [Double] -> [Double] -> [Double]
stochasticAntiAlias rs (lsll:lsl:lsc:lsr:lsrr:[]) = []
stochasticAntiAlias (r:rLst) (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)
    = ( cubInterpolate lsll lsl lsc lsr lsrr lsrrr r )
          : stochasticAntiAlias rLst (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)

rand list adalah daftar variabel acak dalam rentang [0,1].

Meskipun pendekatan Anda secara teoritis benar (dan perlu sedikit dimodifikasi untuk fungsi non-monoton), sangat sulit untuk menghitung kebalikan dari fungsi generik. Seperti yang Anda katakan Anda harus berurusan dengan poin cabang dan pemotongan cabang, yang bisa dilakukan, tetapi Anda benar-benar tidak mau.

Seperti yang telah Anda sebutkan, sampel sampel biasa set poin yang sama dan karenanya sangat rentan terhadap perkiraan yang buruk di daerah di mana sampel tidak (bahkan jika kriteria Nyquist terpenuhi). Dalam hal ini, pengambilan sampel untuk periode yang lebih lama juga tidak membantu.

Secara umum, ketika berhadapan dengan fungsi kerapatan probabilitas dan histogram, itu ide yang jauh lebih baik untuk berpikir dalam hal pengambilan sampel stokastik daripada pengambilan sampel biasa (lihat jawaban terkait untuk pengantar). Dengan mengambil sampel secara stokastik, Anda dapat memastikan bahwa setiap poin memiliki probabilitas yang sama untuk menjadi "hit" dan merupakan cara yang jauh lebih baik untuk memperkirakan pdf.

$f(x)=\sin(20\pi x) + \sin(100\pi x)$$f_s=1000$$f_N=100$$1000$ sampel (distribusi seragam) per detik (saya tidak menggunakan Hz di sini, karena itu menyiratkan arti yang berbeda) selama 30 detik memberikan plot di sebelah kanan (binning yang sama).

Anda dapat dengan mudah melihat bahwa meskipun berisik, ini adalah perkiraan yang jauh lebih baik untuk PDF yang sebenarnya daripada yang di kanan yang menunjukkan nol dalam beberapa interval dan kesalahan besar di beberapa lainnya. Dengan memiliki waktu pengamatan yang lebih lama, Anda dapat menurunkan varians yang ada di sebelah kanan, akhirnya konvergen ke PDF yang tepat (garis hitam putus-putus) dalam batas pengamatan besar.

Estimasi Kepadatan Kernel

Salah satu cara untuk memperkirakan PDF suatu bentuk gelombang adalah dengan menggunakan penduga kepadatan kernel .

$x(n)$$K(\mathrm{x})$$\delta(\mathrm{x} - x(n))$$\hat{P}$

Pembaruan: Informasi tambahan yang menarik.

$x(n)$$n= 0,1,...,N-1$$X(k)$

$X(k)$$e^{\jmath 2\pi n k /N}$

Jadi coba tebak apa yang Anda inginkan adalah untuk menggabungkan semua PDF dari masing-masing komponen Fourier:

$X(k)$$x(n)$

Namun, lebih banyak pemikiran diperlukan!

Seperti yang Anda tunjukkan dalam salah satu komentar Anda, akan sangat menarik untuk dapat menghitung histogram dari sinyal yang direkonstruksi hanya dengan menggunakan sampel dan PDF dari fungsi sinc yang menginterpolasi sinyal bandlimited. Sayangnya, saya tidak berpikir ini mungkin karena histogram sinc tidak memiliki semua informasi yang dimiliki oleh sinyal itu sendiri; semua informasi tentang posisi domain waktu tempat setiap nilai dijumpai hilang. Ini membuat tidak mungkin untuk memodelkan bagaimana versi skala dan waktu-tertunda dari sinc akan dijumlahkan, yang Anda inginkan untuk menghitung histogram dari versi "kontinyu" atau sampel sampel dari sinyal tanpa benar-benar melakukan pengambilan sampel.

Saya pikir Anda dibiarkan dengan interpolasi sebagai pilihan terbaik. Anda memang menunjukkan beberapa masalah yang mencegah Anda dari ingin melakukan ini, yang saya pikir dapat diatasi:

• Biaya komputasi: Ini tentu saja selalu menjadi perhatian relatif, tergantung pada aplikasi spesifik yang Anda inginkan. Berdasarkan tautan yang Anda poskan ke galeri rendering yang telah Anda kumpulkan, saya berasumsi Anda ingin melakukan ini untuk visualisasi sinyal audio. Apakah Anda tertarik pada ini untuk aplikasi real-time atau offline, saya akan mendorong Anda untuk membuat prototipe interpolator yang efisien dan melihat apakah itu benar-benar terlalu mahal. Resampling Polyphase adalah cara yang baik untuk melakukan ini yang fleksibel (Anda dapat menggunakan faktor rasional apa pun).

• $\pi$

Anda perlu merapikan histogram (ini akan menghasilkan hasil yang sama seperti menggunakan metode kernel). Bagaimana tepatnya smoothing harus dilakukan perlu eksperimen. Mungkin juga bisa dilakukan dengan interpolasi. Selain penghalusan, saya yakin Anda juga akan mendapatkan hasil yang lebih baik jika Anda mengganti bentuk gelombang sedemikian rupa sehingga frekuensi pengambilan sampel 'secara signifikan lebih tinggi' daripada frekuensi tertinggi dalam input Anda. Ini akan membantu dalam kasus 'rumit' di mana gelombang sinus terkait dengan frekuensi pengambilan sampel sedemikian rupa sehingga hanya beberapa tempat sampah dalam histogram yang terisi. Jika diambil secara ekstrim, laju sampel yang cukup tinggi akan memberi Anda plot yang bagus tanpa perataan. Jadi upampling dikombinasikan dengan beberapa jenis smoothing akan menghasilkan plot yang lebih baik.

Anda memberikan contoh nada 1kHz, di mana plotnya tidak seperti yang Anda harapkan. Ini proposal saya (kode Matlab / Oktaf)

pixels_vertical = 100;
% This needs to be tuned to your configuration and acceptance
upsampling_factor = 16*(pixels_vertical/100); 
fs_original = 48000;
fsine = 1000; % in Hz
fs_up = upsampling_factor*fs_original;
duration = 1; % in seconds
x = sin(2*pi*fsine*[0:duration*fs_up]/fs_up);
period_in_samples = fs_up/fsine;
hist_points = linspace(-1,1,pixels_vertical);
istart = 1;
iend   = period_in_samples;
pixel_values = hist(x(istart:iend), hist_points);
% smooth pixel values
[b,a] = butter(2,0.2);
pixel_values_smooth = filtfilt(b,a,pixel_values);
figure;hold on;
plot(hist_points, pixel_values);
plot(hist_points, pixel_values_smooth,'r');

Untuk nada 1000Hz Anda, Anda mendapatkan ini

Yang perlu Anda lakukan adalah menyetel ekspresi upsampling_factor sesuai keinginan Anda.

Masih belum 100% tahu persis apa kebutuhan Anda. Tetapi menggunakan prinsip upsampling dan smoothing di atas Anda akan mendapatkan ini untuk nada 1kHz (dibuat dengan Matlab). Perhatikan bahwa dalam histogram mentah ada banyak tempat sampah dengan nol hit.