空間周波数の操作１（ハニング窓）
画像データをフーリエ変換して、空間周波数成分（フーリエスペクトル）を得る方法について、複数回にわたってご紹介したいと思います。
まずはサンプルプログラムを動かすための白黒（グレースケール）の正方形の画像を準備してください。縦と横のサイズは同じで、かつ２のべき乗の画像を準備してください。（フーリエ変換を高速に行うために、ファイルの縦横のサイズは２のべき乗がよいとされていますが、そうでなくても大きな問題はないはずです）
画像の準備ができない場合は、以下のうさぎの画像をお使いください。これはPTBと一緒に配布されている画像でAlphaImageDemoなどで使われています。これをPhotoshopで512x512のサイズにトリミングして、かつグレースケールに変換しています。
さて、フーリエ変換を施す場合には、一般的に画像にハニング窓などを使ってフィルタリングを行い、画像の枠の効果を減少させる処理を行います。MATLABでは、hann関数を使うことでこの処理を行うことができるようなのですが、そのためにはSignal Processing Toolbox が必要で、私はこのToolboxを所有していないため手動で行いました。ただしここでご紹介する方法が適切なのか自信がないところもありますので、各人の責任のもとで動作確認を行っていただければと思います。
ここではハニング窓を枠にそうように四角形でかける場合と、円形にかける場合の２種類の方法をご紹介しています。mesh関数を使って、それぞれのハニング窓の形状を確認することができます。四角形の場合は、図形の対角線上の輝度が若干高くなることにご注意ください。
正規化については、Cole and Wilkins (2013) を参考にしています。ただし、平均輝度を128にしています。（彼らの論文では125です）
サンプルプログラム（spectrum1）を動かすには、上のうさぎの画像をプログラムと同じフォルダに保存の上、実行してください。プログラムを実行すると、ハニング窓をかけた結果が別ウィンドウで表示され、プログラムと同じフォルダ内に、正規化された画像がjpg形式で出力されます。
imagesc関数は、画面上で見やすくなるように画像データ（輝度値）を自動で調整します。imagesc関数で出力されたfigureと、imwriteで画像ファイルとして出力したものとは、見た目が異なりますのでご注意ください。
function spectrum1close all;clear all;imgFileName = 'grayRabbit.jpg'; % 画像のファイル名im1 = imread(imgFileName); % 画像のデータを読み込む。% 画像のファイル名を分解[pathstr, fileName, ext] = fileparts(imgFileName); %% 必要に応じて正規化im2 = double(im1);imMean=mean(im2(:)); % 平均輝度imSD=std(im2(:)); % 輝度の標準偏差im2=(im2-imMean)./imSD; % 平均を０、標準偏差を１に。im2 = im2.*25; % 標準偏差の変更（25は適当）im2 = im2+128; % 平均値の変更（128は256階調の半分）% 正規化しないときは次の行をコメントアウト%im2 = im1;% 正規化した画像をjpgファイルとして出力。imwrite(uint8(im2), [fileName '-norm.jpg'], 'jpeg');%% 座標の設定imSize = size(im1);halfImgSize = imSize(2)/2; % 画像が正方形、つまり、imSize(1) = imSize(2)を前提としている。[X, Y] = meshgrid(-halfImgSize:(halfImgSize-1), halfImgSize:-1:-(halfImgSize-1)); % 行がY方向、列がX方向であることに注意。Y軸は上が正であることに注意。[TH, R] = cart2pol(X,Y); % X,Yを極座標へ変換%% 出力結果を表示するfigureの設定ScreenSize = get(0, 'ScreenSize');figure('Position', [ScreenSize(3)./4, ScreenSize(4)./4, 1000, 800]); % 画像の位置%% 正規化した画像を表示subplot(2,3,1);ShowGrayScale(im2);%% 画像データにハニング窓（四角形）をかける absXY = max(abs(X), abs(Y)); % absXYは0からhalfImgSizeを取り得る。% 画像の中心部では、absXYは0で、hanningの値は1となる。% absXYが最大(halfImgSize）のとき、hanningの値は0となる。hanning1 = 0.5 - 0.5 * cos(2*pi*(halfImgSize/2+absXY./2)/halfImgSize);im3 = double(im2) .* hanning1;subplot(2,3,2);ShowGrayScale(im3);%% 画像データにハニング窓（円形）をかけるR2 = min(R, halfImgSize); % halfImgSizeを超えない半径のデータを作る。% R2は0からhalfImgSizeを取り得る。% 画像の中心部では、R2は0で、hanningの値は1となる。% R2が最大(halfImgSize）のとき、hanningの値は0となる。hanning2 = 0.5 - 0.5 * cos(2*pi*(halfImgSize/2+R2./2)/halfImgSize);im4 = double(im2) .* hanning2;subplot(2,3,3);ShowGrayScale(im4);figure;mesh(hanning1);figure;mesh(hanning2);function ShowGrayScale(imageData)% 画像をグレースケールで表示する。imagesc(imageData);axis('image');axis('off');colormap(gray);