Перейти к содержанию
    

Спектральный анализ Фурье

18 minutes ago, dsl2640free said:

кусок кода, детектирующий гармоники в составе сигнала

Насколько я понял - это обычный ДПФ...

29 minutes ago, ViKo said:

Сюда ходи.

Благодарю. Вроде бы неплохо изложено в интересующем меня месте. Попробую еще раз зайти самостоятельно.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Ну всё, разобрался. Всем спасибо.

Для dsPIC33 одна комплексная "бабочка" требует 20 машинных циклов из которых собственно mac-операции съедают 12, остальное уходит на адресацию исходных операндов и выгрузку результата. Оценочно, для dsPIC33E на 1024-точечный БПФ потребуется примерно 1,6 мс при максимальной скорости его работы (70 MIPS).

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Особенности указателей в DSP модуле dsPIC33 увеличили длину кода "бабочки" до  23 инструкций.

Получилось так:

; W10=pA
; W11=pB
; W8=pW
; A=A+B
; B=(A-B)W
BtflCmplx:
	mov		#0x7FFF, W4				; W4=0.99999
	clr		A, [W10]+=2, W6				; W6=Ar, [W10]->Ai
	clr		B, [W8]+=2, W5, [W11]+=2, W7		; W7=Br, [W11]->Bi, W5=Wr [W8]->Wi
;--- perform A
	mac		W4*W6, A, [W10]-=2, W6			; ACCA=+(1*Ar), W6=Ai, [W10]->Ar
	mac		W4*W7, A, [W11]-=2, W7			; ACCA=+(1*Br), W7=Bi, [W11]->Br
	sac.r		A, #-1, [W13++]				; Ar(new)=ACCA.rnd->buf
	mac		W4*W6, B, [W10]+=2, W6			; ACCB=+(1*Ai), W6=Ar, [W10]->Ai
	mac		W4*W7, B, [W11]+=2, W7			; ACCB=+(1*Bi), W7=Br, [W11]->Bi
	sac.r		B, #-1, [W13++]				; Ai(new)=ACCB.rnd->buf
;--- perform B
	mpy		W5*W6, A, [W10]-=2, W6			; ACCA=Wr*Ar, W6=Ai, [W10]->Ar
	msc		W5*W7, A, [W8]-=2, W5, [W11]-=2, W7	; ACCA=-(Wr*Br), W5=Wi, [W8]->Wr, W7=Bi, [W11]->Br
	msc		W5*W6, A				; ACCA=-(Wi*Ai)
	mac		W5*W7, A, [W8]+=2, W5			; ACCA=+(Wi*Bi), W5=Wr, [W8]->Wi
	sac.r		A, #-1, [W13++]				; Br(new)=ACCA.rnd->buf
	mpy		W5*W6, B, [W10]+=2, W6			; ACCB=Wr*Ai, W6=Ar, [W10]->Ai
	msc		W5*W7, B, [W8], W5, [W11]+=2, W7	; ACCB=-(Wr*Bi), W5=Wi, W7=Br, [W11]->Bi
	mac		W5*W6, B				; ACCB=+(Wi*Ar)
	msc		W5*W7, B				; ACCB=-(Wi*Br)
	sac.r		B, #-1, [W13]				; Bi(new)=ACCB.rnd->buf
	mov		[W13--], [W11--]			; save Bi
	mov		[W13--], [W11]				; save Br
	mov		[W13--], [W10--]			; save Ai
	mov		[W13], [W10]				; save Ar
	return

Итого - "бабочка" выполняется за 330 нс.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Спойлер

 

Ниже приведен пример более сложной процедуры, которая использует DSP особенности dsPIC33E.

Спойлер
#include <dsp.h>;
#define FFT_BLOCK_LENGTH 1024 //Number of frequency points in the FFT
#define LOG2_BLOCK_LENGTH 10 //Number of "Butterfly" Stages in FFT processing should be related to FFT_BLOCK as in 2^9=512
#define AUDIO_FS 44100//Sampling frequency of audio signal captured by the mic
int16_t peakFrequencyBin;
uint16_t ix_MicADCbuff;
uint16_t peakFrequency;
extern const fractcomplex twiddleFactors[FFT_BLOCK_LENGTH/2]__attribute__ ((space(auto_psv), aligned (FFT_BLOCK_LENGTH*2)));
fractcomplex sigCmpx[FFT_BLOCK_LENGTH] __attribute__ ((section (".ydata, data, ymemory"), aligned (FFT_BLOCK_LENGTH * 2 *2)))={{0}};
void readMic(void)//Sample microphone input
{
ADC1_ChannelSelectSet(ADC1_AI_MIC);
ix_MicADCbuff=0;
for(ix_MicADCbuff=0;ix_MicADCbuff&lt;FFT_BLOCK_LENGTH;ix_MicADCbuff++)
{

//Insert delay subroutine here
ADC1_SamplingStop();
while(!ADC1_IsConversionComplete()){}
sigCmpx[ix_MicADCbuff].real = ADC1_Channel0ConversionResultGet();
sigCmpx[ix_MicADCbuff].imag = 0;
}
}

void signalFreq(void)//Detect the dominant frequency of the audio picked by the microphone
{
readMic();
FFTComplexIP (LOG2_BLOCK_LENGTH, &amp;sigCmpx[0], (fractcomplex *) __builtin_psvoffset(&amp;twiddleFactors[0]), (int) __builtin_psvpage(&amp;twiddleFactors[0]));// Perform FFT operation
BitReverseComplex (LOG2_BLOCK_LENGTH, &amp;sigCmpx[0]);// Store output samples in bit-reversed order of their addresses
SquareMagnitudeCplx(FFT_BLOCK_LENGTH, &amp;sigCmpx[0], &amp;sigCmpx[0].real);//Compute the square magnitude of the complex FFT output array so we have a Real output vector
VectorMax(FFT_BLOCK_LENGTH/2, &amp;sigCmpx[0].real, &amp;peakFrequencyBin);//Find the frequency Bin ( = index into the SigCmpx[] array) that has the largest energy
peakFrequency = peakFrequencyBin*(AUDIO_FS/FFT_BLOCK_LENGTH); //Compute the frequency (in Hz) of the largest spectral component
}

 

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Этот код использует не фичи dsPIC, а БИБЛИОТЕКУ, написанную кем то на АСМе. Потому что написать на голом Си работу с DSP ядром этой платформы совершенно невозможно. Полагаю, что и на аналогичных архитектурах DSP других производителей - тоже. Даже не представляю себе, как заставить компилятор управлять указателями и загрузкой по ним в инструкциях DSP...

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Строго говоря, БПФ тут реализовано двумя строками (даже одной):

On 1/27/2020 at 3:50 PM, dsl2640free said:

FFTComplexIP (LOG2_BLOCK_LENGTH, &amp;sigCmpx[0], (fractcomplex *) __builtin_psvoffset(&amp;twiddleFactors[0]), (int) __builtin_psvpage(&amp;twiddleFactors[0]));// Perform FFT operation
BitReverseComplex (LOG2_BLOCK_LENGTH, &amp;sigCmpx[0]);// Store output samples in bit-reversed order of their addresses

сказать, как и насколько функция  FFTComplexIP (...) учитывает особенности dsPIC33, на мой взгляд, нельзя. 

если на асме прилично оптимизировать бабочку(правильные предвыборки и логичные смещения указателей), то дальше ускорить БПФ практически некуда

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

On 1/27/2020 at 5:58 PM, my504 said:

Этот код использует не фичи dsPIC, а БИБЛИОТЕКУ, написанную кем то на АСМе. Потому что написать на голом Си работу с DSP ядром этой платформы совершенно невозможно. Полагаю, что и на аналогичных архитектурах DSP других производителей - тоже. Даже не представляю себе, как заставить компилятор управлять указателями и загрузкой по ним в инструкциях DSP...

 Выделяет характерные паттерны и в путь. Когда писал (немного) для Shark - тоже удивлялся, он и инструкции правильно распараллеливал, и аппаратные циклы использовал. Тоже самое сейчас вижу код под Neon создаётся- править после практически нечего. А скоро думаю AI направят на все это дело и программистам останется только ТЗ грамотно писать.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Присоединяйтесь к обсуждению

Вы можете написать сейчас и зарегистрироваться позже. Если у вас есть аккаунт, авторизуйтесь, чтобы опубликовать от имени своего аккаунта.

Гость
К сожалению, ваш контент содержит запрещённые слова. Пожалуйста, отредактируйте контент, чтобы удалить выделенные ниже слова.
Ответить в этой теме...

×   Вставлено с форматированием.   Вставить как обычный текст

  Разрешено использовать не более 75 эмодзи.

×   Ваша ссылка была автоматически встроена.   Отображать как обычную ссылку

×   Ваш предыдущий контент был восстановлен.   Очистить редактор

×   Вы не можете вставлять изображения напрямую. Загружайте или вставляйте изображения по ссылке.

×
×
  • Создать...