3

FFTと周波数の変更およびforループのベクトル化

あいさつ

以下のコードでfftとフーリエ級数展開FORループの組み合わせを使用して信号の周波数を増減できますが、信号/配列が大きい場合は非常に遅くなります（1x44100の配列は完了するのに約2分かかります） ）forループに関係していることは確かですが、パフォーマンスを向上させるためにベクトル化する方法が正確にはわかりません。これは、3〜6分の長さのオーディオ信号で使用されることに注意してください。1x44100アレイはわずか1秒で、完了するまでに約2分かかります

推奨事項

%create signal
clear all, clc,clf,tic
x= linspace(0,2*pi,44100)';

%Used in exporting to ycalc audio file make sure in sync with above
freq_orig=1;
freq_new=4
vertoff=0;
vertoffConj=0;
vertoffInv=0;
vertoffInvConj=0;
phaseshift=(0)*pi/180 ; %can use mod to limit to 180 degrees

y=sin(freq_orig*(x));
[size_r,size_c]=size(y);

N=size_r; %to test make 50
T=2*pi;
dt=T/N;
t=linspace(0,T-dt,N)';
phase = 0;
f0 = 1/T; % Exactly, one period

y=(y/max(abs(y))*.8)/2; %make the max amplitude here
C = fft(y)/N; % No semicolon to display output


A = real(C);
B = imag(C)*-1; %I needed to multiply by -1 to get the correct sign

% Single-Sided (f >= 0)
An = [A(1); 2*A(2:round(N/2)); A(round(N/2)+1)]; 
Bn = [B(1); 2*B(2:round(N/2)); B(round(N/2)+1)];

pmax=N/2;
ycalc=zeros(N,1); %preallocating space for ycalc
w=0;

for p=2:pmax
               %
       %%1 step) re-create signal using equation
        ycalc=ycalc+An(p)*cos(freq_new*(p-1).*t-phaseshift)
+Bn(p)*sin(freq_new*(p-1).*t-phaseshift)+(vertoff/pmax);
        w=w+(360/(pmax-1)); %used to create phaseshift
        phaseshift=w;

end;

fprintf('\n- Completed in %4.4fsec or %4.4fmins\n',toc,toc/60);

subplot(2,1,1), plot(y),title('Orginal Signal');
subplot(2,1,2),plot(ycalc),title('FFT new signal');

これは、誰かが出力を見たい場合のプロットの写真です。これは正しいです。FORループは本当に本当に遅いです。

ここに画像の説明を入力してください

4

1 に答える 1

1

基本的に周波数領域で信号を上方にシフトしているように見え、「直列展開」はシフトされたバージョンに逆DFTを実装しているだけです。ご覧のとおり、ナイーブなiDFTは非常に遅くなります。そのループ全体をifftの呼び出しに変更してみてください。そうすれば、大幅なスピードアップが得られるはずです。

于 2011-04-02T23:35:48.667 に答える