読者です 読者をやめる 読者になる 読者になる

wavread/wavwriteからaudioread/audiowriteへの移行

Matlab R2012b以降のMatlabではwavread/wavwriteではなくaudioread/audiowriteが推奨されるようになり、エディタ画面でも黄色の警告メッセージ(「deprecatedだよ」)が表示されていました。Matlab R2014bのエディタ画面では赤色のエラーメッセージ(「remov…

Matlab 2014aでドキュメント見られない問題が解決

Matlab 2014aでドキュメントを見ようとすると固まってしまう問題が発生して困っていました。Mavericksにしてからだったかなぁ。helpコマンドは通常通り動くのですが、docコマンドを実行すると、くるくるビーチボール状態になってしまって強制終了を余儀なく…

HomebrewでOctaveのビルドが(ようやく)通った

Mac OS 10.9 MavericksにしてからHomebrewでもMacPortsでもOctaveのビルドが通らず、そのうちにOctave 3.8.0が登場して、それもビルドできず、原因もよく分からない、という状態が続いていた(エラーメッセージを読むのがだるくてサボっていたという面もある…

Matlab/Octaveの判別

Octaveだけに存在するOCTAVE_VERSIONというビルトイン関数があるかどうかを判定の材料に使用します。existは、引数に与えた名称がビルトイン関数であれば5を戻す関数です。(その他、変数なら1、ファイル名なら2、.mexや.octなら3、ディレクトリなら7など、様…

OctaveをMountain Lionにインストールする

OS X 10.8 Mountain LionにHomebrewを使ってOctaveをインストールする方法を、Mountain Lion: Scientific Development SetupとGet Mountain Lion and Homebrew to Be Happyを参考にまとめてみた。自分がつまづいた点もあるので修正点なども加筆した。 Xcode…

モンテカルロ法で円周率の計算

今年度、Octaveで簡単な統計計算をする授業をやっているのですが、その中で大数の法則と中心極限定理の説明をするための一例として、モンテカルロ法で円周率を計算する練習問題を出題しました。座標軸上(0,0)-(1,1)に1辺が長さ1の正方形があり、それに接する…

Octave 3.4.0バイナリ・パッケージをLionで動かす

2月13日(月)に更新されたREADMEファイルの抄訳です。Mac OS X Lion 10.7.3でOctave-Forgeのバイナリに付属しているGnuplotを動かそうとすると以下のようなメッセージが出てプロットがされません。 dyld: Library not loaded: /usr/X11/lib/libfreetype.6.dyl…

fftfiltを使うときに注意すること

Matlab 2011bで、単純な低域フィルタ [0.5 1 0.5] とインパルス的な [5 1 1 1 1] をfftfiltを使って畳み込みしてみたところ、fftfiltをそのまま使った場合、尻尾が切られてしまいます。 >> fftfilt([0.5 1 0.5], [5 1 1 1 1]) ans = 2.5000 5.5000 4.0000 2.…

数字を読み上げる音声ファイルを作る

主観評価実験をするときに、刺激を逐次呈示して、それに対する評価を記述してもらう方法があります。10個程度の刺激ならあまり問題にはならないのですが、刺激が多くなると、いま何番目の刺激なのか、いま使っている回答欄は正しいのか、など不安になること…

複数ファイルを逐次処理する

似た種類のファイルがたくさんあって、それらをひとつずつ読み込んで処理したい、ということがときどきあります。たとえば、複数のWAVファイルすべてをノーマライズしたいとか、たくさんあるテキストファイルの1行目だけを読み込んで別ファイルに書き出した…

音響数値解析のお勉強(有限要素法)

若手の音響研究者たちが作成しているOpenAcousticsというドキュメント+ソースコード群があります。このサイトでは音響数値解析のソースコードだけでなく、充実した日本語の入門書が無料公開されています。だいぶ前に教えてもらっていたのですが、重い腰を上…

MacPorts版OctaveをLionで動かすために

Lion導入に合わせて/usr/local以下をすべて削除し、LaTeXなどコマンドラインから使用するソフトの多くにMacPortsを利用しはじめています。数値計算ソフトのOctaveはこれまでバイナリ版を利用してきましたが、Octave 3.2.3はMac OS X 10.6.5から(比較的簡単…

一様乱数で正規乱数をつくる

先日、学生から「Octaveではrandn()で正規分布に従う乱数が作れるが、C/C++にはそのようなものが用意されていないようで困っている」との相談がありました。聞いてみると、厳密に正規分布に従う必要はなく、ゼロ近傍の発生確率が高く、ゼロから離れるに従っ…

Project Euler No.14

数日ぶりにProject Euler。第14問。 The following iterative sequence is defined for the set of positive integers:n -> n/2 (n is even) n -> 3n + 1 (n is odd)Using the rule above and starting with 13, we generate the following sequence:13 -> 4…

Octave 3.4.0 for Mac OS Xのgnuplot同梱版

SourceForgeで公開されているMac OS X用のOctave 3.4.0の2011-04-08版にはgnuplotが入っておらず、plotなどが使えません。それに応える形で、現在のメンテナーであるJulien Salort氏のサイトでgnuplot同梱版が公開されました。2011-04-21版が本日時点での最…

配列を宣言してから中身を入れると速い

「R in a Nutshell」を読み進んでいます。第8章に、配列というのはあらかじめメモリを確保してから中身を入れた方が速いということが説明されています。たとえば create.vector.of.ones <- function(n) { return.vector <- NA; for (i in 1:n) { return.vect…

Project Euler No.13

Project Euler第13問。 Work out the first ten digits of the sum of the following one-hundred 50-digit numbers.(数表省略)以下に100桁の数が50個ある。すべて加算したときに、最上位の10桁の数は何か。 これはMatlabではラクチン。この程度の桁数なら問…

Project Euler No.11

Project Eulerの11問目。 In the 2020 grid below, four numbers along a diagonal line have been marked in red.(数表省略)The product of these numbers is 26 × 63 × 78 × 14 = 1788696.What is the greatest product of four adjacent numbers in any d…

Project Euler No.10

さてさて第10問目のProject Euler。 The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.Find the sum of all the primes below two million.10未満の素数の和は17である。2百万未満の全素数の和はいくつか。 Matlabにはprimes関数があるので楽勝です。…

Project Euler No.9

次々行くよ。Project Eulerの9問目。 A Pythagorean triplet is a set of three natural numbers, a b c, for which,a2 + b2 = c2 For example, 32 + 42 = 9 + 16 = 25 = 52.There exists exactly one Pythagorean triplet for which a + b + c = 1000. Find…

Project Euler No.8

引き続きProject Eulerの8問目。 Find the greatest product of five consecutive digits in the 1000-digit number.73167176531330624919225119674426574742355349194934 96983520312774506326239578318016984801869478851843 8586156078911294949545950173…

Project Euler No.7

MatlabでProject Eulerを解く第7問。 By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.What is the 10001st prime number?最初の6つの素数を並べると2、3、5、7、11、13となる。6番目の素数は13だと…

Project Euler No.5

引き続きProject EulerをMatlab/Octaveで解いてます。 2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.What is the smallest positive number that is evenly divisible by all of the number…

Project Euler No.4

Project Eulerの第4問は「数字の回文」。 A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 × 99.Find the largest palindrome made from the product of two 3-digit nu…

Project Euler No.3

Project Eulerの問題をMatlab/Octaveで解いてみる3問目。 The prime factors of 13195 are 5, 7, 13 and 29. What is the largest prime factor of the number 600851475143 ?13195の素因数は5、7、13、29である。600851475143の素因数のうち最大のものは何…

Project Euler No.2

Project Eulerの問題をMatlabで解いてみる2問目。 Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...By considering…

Project Euler No.1

Project Eulerはコンピュータで解くための問題集を公開しているサイトです。その第1問目をMatlabでやってみました。Octaveでも動きます。 If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of the…

相互相関と畳み込み

相互相関と畳み込みは処理が似ているので、どちらかを使ってもう一方を計算できないかと考えた。で、Matlabで実験したところ、xcorr(a, b)はconv(a, flipud(b))と同じだということがわかった。convの代わりにFFTを用いて処理してみたが*1、xcorrのほうが若干…

Snow Leopard(およびLeopard)でOctave.appを動かす

Octave Forgeで配布されているOctaveのMac OS Xバイナリがありますが、2009年10月に作成されたものが現時点での最新版で、Mac OS X 10.6.5 Snow Leopardおよび10.5.8に対応していません。そのため、Gnuplotを用いたグラフのプロットなどができなくなっていま…

Octaveによる主成分分析の計算

Matlabにはprincompという主成分分析を計算する関数があるけれど、Octaveには見あたらなかったので、英語版Wikipediaに書かれていたやり方どおりに作ってみました。 [pc, score, latent] = princomp(A); みたいに使えます。pcが主成分の係数、scoreが主成分…

Matlab/Octaveで平均律

単純なプログラムですが、MIDIノート番号を周波数に、周波数をMIDIノート番号とピッチベンド値(およびセント値)に変換するものを置いておきます。特に後者はあまり見かけないと思うので。デフォルト設定は、基準ピッチが440Hz、ベンドレンジは±2半音です。…

遺伝的アルゴリズム的な何かでフィルターを設計してみる(でもあまりうまくいかない)

アストロノーカで遊んで興味を持ったのが遺伝的アルゴリズム。その考え方を使ってフィルタ設計でもできないかと思ったので、とりあえずやってみました。中に出てくるbiquad_*は、Robert Bristow-Johnsonによるフィルター設計レシピの実装で、中心周波数やゲ…

ラウドネスの近似値計算

コンピュータの中にあるWAVファイルどうしの音量比較をしたいときはけっこうあります。音量の指標として「ラウドネス(sone)」があり、以前その計算プログラムも書いたのですが、その音を鳴らしたときの音圧値(dBSPL値)が分からないと正確な値が得られな…

Matlabで2次元プロット

論文に載せる2次元プロットを作成するとき、僕はMatlabを使って以下のような感じで書いています。多次元尺度法の分析結果を図にするときに相性がいいですよ。背景になる部品から順番に重ねていく必要があるので、まずは縦軸と横軸を書きます。grid onでグリ…

Matlab R2009b

Matlab R2009bからMac OS X版が64ビット対応になっていて良い感じ。もしかしたらそれ以前から64ビット対応していたのかもしれないけど、Snow Leopardで分かるようになりました。計算速度が大きく変化したような気はしませんが、64ビットだと言うだけで気持ち…

音の周波数特性を描くMatlabプログラム

音の周波数特性を描くMatlabプログラムを改善したので、メモ的に置いておきます。今回の目玉は「1/3オクターブごとの平均値を描くようにした」です。500-2000Hz帯の平均値が0dBになるようにスケーリングしています。

遺伝的アルゴリズムの簡単なテスト

以前から遺伝的アルゴリズムというものに興味を持っていたのですが、本を読むだけじゃなく実際に手を動かさなければ、ということでプログラムを書いてみました。勉強中のLISPで書きたいところでしたが、いちばん慣れているMatlabを使ってみました。の最大値…

シェッフェの一対比較法(中屋変法)のMatlabコード

研究で必要になったので、佐藤信『統計的官能検査法』(日科技連 1985年)を参考に、シェッフェの一対比較法(中屋変法)による分散分析を行うMatlabスクリプトを作りました。実のところ「参考にした」どころではなく、書いてあった手計算のアルゴリズムをそ…

Matlabソースの高速化のために

Matlabソースの高速化のためにできることのまず第一は、ループ(forやwhile)を少なくし、行列内の各々の要素にアクセスする頻度を少なくすることです。そのためにはコードのベクトル化(code vectorization)をしてやる必要があります。実はこの手法(code ve…

ロジットとプロビット

ロジットしか使ったことなかったのですが、プロビットも同じ目的で使えるのですね。どちらも独立変数が連続値で従属変数がブール値のときに、シグモイド関数を重回帰分析ではめる分析手法です。SPSSのページによると、違いは「ロジスティック関数を使うか正…

インパルス応答の分析(残響時間)

ほっしくんに指摘されたので、間違った情報(10月13日のエントリー後半)を修正すべく、プログラムを書きました。残響時間(RT60)は、「信号のエネルギーが60dB減衰するまでの時間」と定義されています。ISOではその定義を発展させて?現実の測定に即した形…

インパルス応答の分析(直間比)

昨日の続き。直接音と間接音のレベル比を計算します。ISO 3382という規格にclarityという物理指標が出てきます。インパルス応答の前半部分(直接音に対応)と後半部分(間接音に対応)のレベル比をdBで求めるもので、例えば以下のような式で表されます。この…

インパルス応答の分析(残響時間を求める)

Matlabで部屋のインパルス応答の分析をしています。分析するのには多種多様な方法があるのですが、まずは残響時間の計算ができて、直接音と間接音の比が求まれば、とりあえずの分析材料にはなります。残響時間を求めるために、以前作成した残響曲線の計算を…

Zoom H2で音源方向認識

Zoom H2というのは4chのPCM録音ができるポータブル・レコーダーなのですが、せっかく4つのマイク位置が固定になっているので、それを使って遊んでみようと思い、Matlabで非常に単純な方法で音源方向の認識を行うMatlabルーチンを作ってみました。H2には前方…

MatlabでMP3を読み書きする

MatlabでMP3を読み書きするルーチンが欲しくなって、作りました。無圧縮ファイルを使うほど音質重視じゃないけど、処理したいファイル数はたくさんある、という時にちょっと便利かもしれません。自分の環境のことしか考えてなかったので、Unix系でLAMEが入っ…