Величина и информация о фазе БПФ

Представление частотного диапазона сигнала несет информацию о величине сигнала и фазе на каждой частоте. Поэтому выход расчета БПФ является комплексным. Комплексное число, имеет действительную часть и мнимую часть, такую что. Величина вычисляется как, и фаза вычисляется как. Можно использовать abs функций MATLAB и угол, чтобы соответственно получить величину и фазу любого комплексного числа.

Используйте аудио пример, чтобы разработать некоторое понимание на том, какую информацию несут величина и фаза сигнала. Для этого загрузите звуковой файл, содержащий 15 секунд музыки акустической гитары. Частота дискретизации звукового сигнала составляет 44,1 кГц.

Fs = 44100;
y = audioread('guitartune.wav');

Используйте fft, чтобы наблюдать содержимое частоты сигнала.

NFFT = length(y);
Y = fft(y,NFFT);
F = ((0:1/NFFT:1-1/NFFT)*Fs).';

Выход БПФ является комплексным вектором, содержащим информацию о содержимом частоты сигнала. Величина говорит вам силу частотных составляющих относительно других компонентов. Фаза говорит вам, как все частотные составляющие выравниваются вовремя.

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

magnitudeY = abs(Y);        % Magnitude of the FFT
phaseY = unwrap(angle(Y));  % Phase of the FFT

helperFrequencyAnalysisPlot1(F,magnitudeY,phaseY,NFFT)

Можно применять обратное преобразование Фурье к вектору частотного диапазона, Y, чтобы восстановить сигнал времени. 'Симметричный' флаг говорит ifft, что вы имеете дело с сигналом времени с действительным знаком, таким образом, это обнулит маленькие мнимые компоненты, которые появляются на обратном преобразовании из-за числовых погрешностей в расчетах. Заметьте, что исходный сигнал времени, y, и восстановленный сигнал, y1, является практически тем же самым (норма их различия находится на порядке 1e-14). Очень небольшая разница между этими двумя происходит также из-за числовых упомянутых выше погрешностей. Вопроизведите и послушайте непреобразованный сигнал y1.

y1 = ifft(Y,NFFT,'symmetric');
norm(y-y1)

ans =

   3.9466e-14

hplayer = audioplayer(y1, Fs);
play(hplayer);

Чтобы видеть эффекты изменения ответа величины сигнала, удалите частотные составляющие выше 1 кГц непосредственно от FFT выход (путем создания величин равными нулю) и слушайте влияние, которое это оказывает на звук звукового файла. Удаление высокочастотных компонентов сигнала упоминается как фильтрация lowpass.

Ylp = Y;
Ylp(F>=1000 & F<=Fs-1000) = 0;

helperFrequencyAnalysisPlot1(F,abs(Ylp),unwrap(angle(Ylp)),NFFT,...
  'Frequency components above 1 kHz have been zeroed')

Верните отфильтрованный сигнал в область времени, использующую ifft.

ylp = ifft(Ylp,'symmetric');

Проигрывайте сигнал. Можно все еще услышать мелодию, но она походит, покрыли ли вы уши (вы фильтруете высокочастотные звуки, когда вы делаете это). Даже при том, что гитары производят примечания, которые являются между 400 и 1 кГц, когда вы проигрываете примечание по строке, строка также вибрирует во множителях основной частоты. Эти более высокие частотные составляющие, называемые гармониками, то, что дает гитаре его конкретный тон. Когда вы удаляете их, вы заставляете звук казаться "непрозрачным".

hplayer = audioplayer(ylp, Fs);
play(hplayer);

Фаза сигнала имеет важную информацию о том, когда вовремя примечания песни появляются. Чтобы проиллюстрировать важность фазы на звуковом сигнале, удалите информацию о фазе полностью путем взятия величины каждой частотной составляющей. Обратите внимание на то, что путем выполнения этого вы сохраняете ответ величины неизменным.

% Take the magnitude of each FFT component of the signal
Yzp = abs(Y);
helperFrequencyAnalysisPlot1(F,abs(Yzp),unwrap(angle(Yzp)),NFFT,[],...
  'Phase has been set to zero')

Верните сигнал во временном интервале и проигрывайте аудио. Вы не можете распознать исходный звук вообще. Ответ величины является тем же самым, никакие частотные составляющие не были удалены на этот раз, но порядок примечаний исчез полностью. Сигнал теперь состоит из группы синусоид все выровненные в равное нулю время. В общем случае искажения фазы, вызванные фильтрацией, могут повредить сигнал на грани рендеринга его неузнаваемый.

yzp = ifft(Yzp,'symmetric');
hplayer = audioplayer(yzp, Fs);
play(hplayer);

Нахождение периодичностей сигнала

Представление частотного диапазона сигнала позволяет вам наблюдать несколько характеристик сигнала, которые или не легко видеть, или не видимые вообще, когда вы смотрите на сигнал во временном интервале. Например, анализ частотного диапазона становится полезным, когда вы ищете циклическое поведение сигнала.

Анализ циклического поведения температуры в офисном здании

Рассмотрите набор измерений температуры в офисном здании в течение зимнего сезона. Измерения проводились каждые 30 минут в течение приблизительно 16,5 недель. Посмотрите в то время доменные данные с осью времени, масштабируемой к неделям. На этих данных могло быть периодическое поведение?

load officetemp.mat
Fs = 1/(60*30);   % Sample rate is 1 sample every 30 minutes
t = (0:length(temp)-1)/Fs;

helperFrequencyAnalysisPlot2(t/(60*60*24*7),temp,...
  'Time in weeks','Temperature (Fahrenheit)')

Почти невозможно знать, существует ли циклическое поведение на офисных температурах путем рассмотрения сигнала временной области. Однако циклическое поведение температуры становится очевидным, если мы смотрим на ее представление частотного диапазона.

Получите представление частотного диапазона сигнала. Если вы строите величину FFT выход с осью частоты, масштабируемой к циклам/неделя, вы видите, что существует две спектральных линии, которые ясно больше, чем какая-либо другая частотная составляющая. Одна спектральная линия находится в 1 цикле/неделя, другой находится в 7 циклах/неделя. Это целесообразно, учитывая, что данные прибывают из терморегулируемого здания на 7-дневном календаре. Первая спектральная линия указывает, что создающие температуры следуют за еженедельным циклом с более низкими температурами в выходные и более высокими температурами в течение недели. Вторая линия указывает, что существует также ежедневный цикл с более низкими температурами во время ночных и более высоких температур в течение дня.

NFFT = length(temp);              % Number of FFT points
F = (0 : 1/NFFT : 1/2-1/NFFT)*Fs; % Frequency vector

TEMP = fft(temp,NFFT);
TEMP(1) = 0; % remove the DC component for better visualization

helperFrequencyAnalysisPlot2(F*60*60*24*7,abs(TEMP(1:NFFT/2)),...
  'Frequency (cycles/week)','Magnitude',[],[],[0 10])

Измерение степени

Функция периодограммы вычисляет БПФ сигнала и нормирует выход, чтобы получить степень спектральная плотность, PSD или спектр мощности, от которого можно измерить степень. PSD описывает, как степень сигнала времени распределяется с частотой, это имеет модули Ватт/Гц. Вы вычисляете спектр мощности путем интеграции каждой точки PSD на интервале частоты, в котором та точка задана (т.е. по пропускной способности разрешения PSD). Модули спектра мощности являются ваттами. Можно считать значения степени непосредственно из спектра мощности, не имея необходимость объединяться на интервале. Обратите внимание на то, что PSD и спектр мощности действительны, таким образом, они не содержат информации о фазе.

Измерение гармоник при Выходе нелинейного усилителя мощности

Загрузите данные, измеренные при выходе усилителя мощности, который имеет третье искажение порядка формы, где выходное напряжение и является входным напряжением. Данные были собраны с частотой дискретизации 3,6 кГц. Вход состоит из синусоиды на 60 Гц с амплитудой единицы. По причине нелинейного искажения необходимо ожидать, что выходной сигнал усилителя будет содержать компонент DC, компонент на 60 Гц, и вторые и третьи гармоники на уровне 120 и 180 Гц.

Загрузите 3 600 выборок усилителя выход, вычислите спектр мощности и постройте результат в логарифмическом масштабе (децибелы-ватты или dBW).

load ampoutput1.mat
Fs = 3600;
NFFT = length(y);

% Power spectrum is computed when you pass a 'power' flag input
[P,F] = periodogram(y,[],NFFT,Fs,'power');

helperFrequencyAnalysisPlot2(F,10*log10(P),'Frequency in Hz',...
  'Power spectrum (dBW)',[],[],[-0.5 200])

График спектра мощности показывает три из четырех ожидаемых peaks в DC, 60, и 120 Гц. Это также показывает еще несколько побочных peaks, которые должны быть вызваны шумом в сигнале. Обратите внимание на то, что гармоника на 180 Гц полностью прокладывается под землей в шуме.

Измерьте степень видимого ожидаемого peaks:

PdBW = 10*log10(P);
power_at_DC_dBW = PdBW(F==0)   % dBW

[peakPowers_dBW, peakFreqIdx] = findpeaks(PdBW,'minpeakheight',-11);
peakFreqs_Hz = F(peakFreqIdx)
peakPowers_dBW

power_at_DC_dBW =

   -7.8873


peakFreqs_Hz =

    60
   120


peakPowers_dBW =

   -0.3175
  -10.2547

Улучшение измерений мощности для сигналов с шумом

Как замечено на графике выше, периодограмма показывает несколько peaks частоты, которые не связаны с сигналом интереса. Спектр выглядит очень шумным. Причина этого состоит в том, что вы только анализировали одну короткую реализацию сигнала с шумом. Повторение эксперимента несколько раз и усреднения удалило бы побочный спектральный peaks и дало бы к более точным измерениям мощности. Можно достигнуть этого усреднения с помощью функции pwelch. Эта функция возьмет большой вектор данных, повредит его в меньшие сегменты заданной длины, вычислит столько же периодограмм, сколько существуют сегменты и составляют в среднем их. Как количество доступных увеличений сегментов, функция pwelch даст к более сглаженному спектру мощности (меньше отклонения) со значениями степени ближе к ожидаемым значениям.

Загрузите большее наблюдение, состоящее из 500e3 точки усилителя выход. Сохраните число точек используемым, чтобы выполнить БПФ как 3 600 так, чтобы пол (500e3/3600) = 138 БПФ был усреднен, чтобы получить спектр мощности.

load ampoutput2.mat
SegmentLength = NFFT;

% Power spectrum is computed when you pass a 'power' flag input
[P,F] = pwelch(y,ones(SegmentLength,1),0,NFFT,Fs,'power');

helperFrequencyAnalysisPlot2(F,10*log10(P),'Frequency in Hz',...
  'Power spectrum (dBW)',[],[],[-0.5 200])

Как замечено на графике, pwelch эффективно удаляет весь побочный peaks частоты, вызванный шумом. Спектральный компонент на уровне 180 Гц, который содержался шум, теперь отображается. Усреднение удаляет отклонение из спектра, и это эффективно дает к более точным измерениям мощности.

Измерение общей средней степени и степени по диапазону частот

Измерение общей средней степени сигнала временной области является легкой и общей задачей. Для выходного сигнала усилителя, y, общая средняя степень вычисляется во временном интервале как:

pwr = sum(y.^2)/length(y) % in watts

pwr =

    8.1697

В частотном диапазоне общая средняя степень вычисляется как сумма степени всех частотных составляющих сигнала. Значение pwr1 состоит из суммы всех частотных составляющих, доступных в спектре мощности сигнала. Значение соглашается со значением мощности, вычисленной выше использования сигнала области времени:

pwr1 = sum(P) % in watts

pwr1 =

    8.1698

Но что, если вы хотели измерить общую степень, доступную по полосе частот? Можно использовать функцию bandpower, чтобы вычислить степень по любому желаемому диапазону частот. Можно передать сигнал временной области непосредственно как вход к этой функции, чтобы получить степень по заданной полосе. В этом случае функция оценит спектр мощности с методом периодограммы.

Вычислите степень по 50 Гц полосе на 70 Гц. Результат будет включать степень на 60 Гц плюс шумовая степень по полосе интереса:

pwr_band = bandpower(y,Fs,[50 70]);
pwr_band_dBW = 10*log10(pwr_band) % dBW

pwr_band_dBW =

    0.0341

Если вы хотите управлять расчетом спектра мощности, используемого, чтобы измерить степень в полосе, можно передать вектор PSD функции bandpower. Например, можно использовать функцию pwelch, как вы сделали прежде, чтобы вычислить PSD и гарантировать усреднение шумовых эффектов:

% Power spectral density is computed when you specify the 'psd' option
[PSD,F]  = pwelch(y,ones(SegmentLength,1),0,NFFT,Fs,'psd');
pwr_band1 = bandpower(PSD,F,[50 70],'psd');
pwr_band_dBW1 = 10*log10(pwr_band1) % dBW

pwr_band_dBW1 =

    0.0798

Нахождение спектральных компонентов

Сила сигнала состоять из одной или нескольких частотных составляющих. Способность наблюдать все спектральные компоненты зависит от разрешения частоты вашего анализа. Разрешение частоты или пропускная способность разрешения спектра мощности заданы как R = Fs/N, где N является продолжительностью наблюдения сигнала. Только спектральные компоненты, разделенные частотой, больше, чем разрешение частоты, будут разрешены.

Анализ системы управления вибрации землетрясения создания

Системы управления Активного массового драйвера (AMD) используются, чтобы уменьшать вибрацию в создании под землетрясением. Активный массовый драйвер помещается в верхний этаж создания и, на основе смещения и ускоряющих измерений этажей создания, система управления отправляет сигналы в драйвер так, чтобы масса переместилась, чтобы ослабить наземные воздействия. Ускоряющие измерения были зарегистрированы на первом этаже трех тестовых структур истории при условиях землетрясения. Измерения были проведены без активной массовой системы управления драйвера (условие разомкнутого цикла), и с активной системой управления (условие замкнутого цикла).

Загрузите ускоряющие данные и вычислите спектр мощности для ускорения первого этажа. Длина векторов данных 10e3, и частота дискретизации составляет 1 кГц. Используйте pwelch с сегментами длины 64 точки данных, чтобы получить пол (10e3/64) = 156 средних значений БПФ и пропускная способность разрешения Фс/64 = 15,625 Гц. Как был показан прежде, усреднение уменьшает шумовые эффекты и дает к более точным измерениям мощности. Используйте 512 точек БПФ. Используя NFFT> N эффективно интерполирует точки частоты, представляющие более подробный график спектра (это достигается путем добавления нулей NFFT-N в конце сигнала времени и взятия БПФ NFFT-точки нуля дополненный вектор).

Разомкнутый цикл и близкие ускоряющие спектры мощности цикла показывают, что, когда система управления активна, ускоряющий спектр мощности уменьшается между 4 и 11 дБ. Максимальное затухание происходит на уровне приблизительно 23,44 кГц. Сокращение на 11 дБ означает, что мощность вибрации уменьшена фактором 12,6. Общая мощность уменьшена от 0,1670 до 0,059 ватт, фактор 2,83.

load quakevibration.mat

Fs = 1e3;                 % sample rate
NFFT = 512;               % number of FFT points
segmentLength = 64;       % segment length

% open loop acceleration power spectrum
[P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,'power');

% closed loop acceleration power spectrum
P1_CL     = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,'power');

helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),...
  'Frequency in Hz','Acceleration Power Spectrum in dB',...
  'Resolution bandwidth = 15.625 Hz',{'Open loop', 'Closed loop'},[0 100])

Вы анализируете данные о вибрации, и вы знаете, что колебания имеют циклическое поведение. Затем как получается, что графики спектра, показанные выше, не содержат резких спектральных линий, типичных для циклического поведения? Возможно, вы пропускаете те линии, потому что они не разрешимы с разрешением, полученным с 64 длинами сегмента точки? Увеличьте разрешение частоты, чтобы видеть, существуют ли спектральные линии, которые не были разрешимы прежде. Сделайте это путем увеличения длины сегмента данных, используемой в pwelch, функционирует к 512 точкам. Это дает к новому разрешению Фс/512 = 1,9531 Гц. В этом случае количество средних значений БПФ сокращено на пол (10e3/512) = 19. Безусловно, существует компромисс между количеством средних значений и разрешением частоты при использовании pwelch. Сохраните количество точек БПФ равным 512.

NFFT = 512;            % number of FFT points
segmentLength = 512;   % segment length

[P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,'power');
P1_CL     = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,'power');

helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),...
  'Frequency in Hz','Acceleration Power Spectrum in dB',...
  'Resolution bandwidth = 1.95 Hz',{'Open loop', 'Closed loop'},[0 100])

Заметьте, как увеличение разрешения частоты позволяет вам наблюдать три peaks относительно спектра разомкнутого цикла и два на близком спектре цикла. Этот peaks не был разрешим прежде. Разделение между peaks на спектре разомкнутого цикла составляет приблизительно 11 Гц, который меньше, чем разрешение частоты, полученное с сегментами длины 64, но больше, чем разрешение, полученное с сегментами длины 512. Циклическое поведение колебаний теперь отображается. Основная частота вибрации на уровне 5,86 Гц, и равномерно расположенный peaks частоты предполагает, что они гармонично связаны. В то время как было уже замечено, что система управления уменьшает полную мощность колебаний, более высокие спектры разрешения показывает, что другой эффект системы управления состоит в том, чтобы отметить гармонический компонент на уровне 17,58 Гц. Таким образом, система управления не только уменьшает вибрацию, но также и приближает ее к синусоиде.

Важно отметить, что разрешение частоты определяется количеством точек сигнала, не количеством точек БПФ. Увеличение числа точек БПФ интерполирует данные о частоте, чтобы предоставить вам больше подробной информации о спектре, но это не улучшает разрешение.

Заключения

В этом примере вы изучили, как выполнить анализ частотного диапазона сигнала с помощью fft, ifft, периодограммы, pwelch, и функций bandpower. Вы изучили сложный характер БПФ и что информация содержится в величине и фазе спектра частоты. Вы видели преимущества использования данных о частотном диапазоне при анализе периодичности сигнала. Вы изучили, как вычисляют общую степень или степень по конкретной полосе частот сигнала с шумом. Вы изучили, как увеличение разрешения частоты спектра позволяет вам наблюдать близко расположенные частотные составляющие, и вы узнали о компромиссе между разрешением частоты и спектральным усреднением.

Дополнительные материалы для чтения

Для получения дополнительной информации о частотном диапазоне анализ см. Signal Processing Toolbox.

Ссылка: Дж.Г. Проукис и Д. Г. Мэнолакис, "цифровая обработка сигналов. Принципы, алгоритмы и приложения", Prentice Hall, 1996.

Приложение

Следующие функции помощника используются в этом примере.