MATLAB FFT xaxis giới hạn rối tung lên và fftshift

Question

Đây là lần đầu tiên tôi sử dụng hàm fft và tôi đang cố gắng vẽ phổ tần số của hàm cosin đơn giản:

f = cos (2 * pi * 300 * t)

Tỷ lệ lấy mẫu là 220500. Tôi đang vẽ một giây của hàm f.

Đây là nỗ lực của tôi:

time = 1; 
freq = 220500; 
t = 0 : 1/freq : 1 - 1/freq; 
N = length(t); 
df = freq/(N*time); 

F = fftshift(fft(cos(2*pi*300*t))/N); 
faxis = -N/2/time : df : (N/2-1)/time; 

plot(faxis, real(F)); 
grid on; 
xlim([-500, 500]); 

Tại sao tôi nhận kết quả kỳ lạ khi tôi tăng tần suất đến 900Hz? Các kết quả lẻ này có thể được cố định bằng cách tăng giới hạn trục x từ 500Hz đến 1000Hz. Ngoài ra, đây có phải là cách tiếp cận chính xác không? Tôi nhận thấy nhiều người khác đã không sử dụng fftshift(X) (nhưng tôi nghĩ rằng họ chỉ làm một phân tích phổ mặt đơn).

Cảm ơn bạn.

Answer 1

Đây là câu trả lời của tôi như đã hứa.

Câu hỏi đầu tiên hoặc liên quan đến lý do bạn "nhận kết quả lẻ khi bạn tăng tần suất lên 900 Hz" có liên quan đến chức năng rescaling cốt truyện của Matlab như được mô tả bởi @Castilho. Khi bạn thay đổi phạm vi của trục x, Matlab sẽ cố gắng hữu ích và rescale trục y. Nếu các đỉnh nằm ngoài phạm vi chỉ định của bạn, MATLAB sẽ phóng to các lỗi số nhỏ được tạo ra trong quá trình. Bạn có thể khắc phục điều này bằng lệnh 'ylim' nếu nó làm phiền bạn.

Tuy nhiên, câu hỏi thứ hai, cởi mở hơn của bạn "đây có phải là cách tiếp cận chính xác không?" yêu cầu thảo luận sâu hơn. Cho phép tôi nói với bạn làm thế nào tôi sẽ đi về việc thực hiện một giải pháp linh hoạt hơn để đạt được mục tiêu của bạn âm mưu một làn sóng cosin.

Bạn bắt đầu với những điều sau:

time = 1; 
freq = 220500; 

Điều này đặt ra một báo động trong đầu tôi ngay lập tức. Nhìn vào phần còn lại của bài đăng, bạn dường như quan tâm đến tần số trong dải phụ kHz. Nếu đúng như vậy, thì tốc độ lấy mẫu này quá mức khi giới hạn Nyquist (sr/2) cho tốc độ này là trên 100 kHz.Tôi đoán bạn có nghĩa là để sử dụng tỷ lệ lấy mẫu âm thanh phổ biến của 22050 Hz (nhưng tôi có thể sai ở đây)?

Dù bằng cách nào, phân tích của bạn hoạt động tốt về số lượng cuối cùng. Tuy nhiên, bạn không giúp mình hiểu cách FFT có thể được sử dụng hiệu quả nhất để phân tích trong các tình huống thực tế.

Cho phép tôi đăng cách tôi thực hiện việc này. Kịch bản sau thực hiện gần như chính xác những gì kịch bản của bạn làm, nhưng mở ra một số tiềm năng mà chúng ta có thể xây dựng. .

%// These are the user parameters 
durT = 1; 
fs = 22050; 
NFFT = durT*fs; 
sigFreq = 300; 

%//Calculate time axis 
dt = 1/fs; 
tAxis = 0:dt:(durT-dt); 

%//Calculate frequency axis 
df = fs/NFFT; 
fAxis = 0:df:(fs-df); 

%//Calculate time domain signal and convert to frequency domain 
x = cos( 2*pi*sigFreq*tAxis ); 
F = abs( fft(x, NFFT)/NFFT ); 

subplot(2,1,1); 
plot( fAxis, 2*F ) 
xlim([0 2*sigFreq]) 
title('single sided spectrum') 

subplot(2,1,2); 
plot( fAxis-fs/2, fftshift(F) ) 
xlim([-2*sigFreq 2*sigFreq]) 
title('whole fft-shifted spectrum') 

Bạn tính trục thời gian và tính số điểm FFT của bạn từ độ dài của trục thời gian. Điều này rất kỳ quặc. Vấn đề với cách tiếp cận này, là độ phân giải tần số của các thay đổi fft khi bạn thay đổi thời lượng tín hiệu đầu vào của bạn, bởi vì N phụ thuộc vào biến "thời gian" của bạn. Lệnh MATLAB MATLAB sẽ sử dụng kích thước FFT khớp với kích thước của tín hiệu đầu vào.

Trong ví dụ của tôi, tôi tính trục tần số trực tiếp từ NFFT. Điều này có phần không liên quan trong bối cảnh của ví dụ trên, khi tôi đặt NFFT bằng số lượng mẫu trong tín hiệu. Tuy nhiên, việc sử dụng định dạng này giúp làm sáng tỏ suy nghĩ của bạn và nó trở nên rất quan trọng trong ví dụ tiếp theo của tôi.

** CHÚ Ý: Bạn sử dụng thực (F) trong ví dụ của mình. Trừ khi bạn có một lý do rất tốt để chỉ trích xuất phần thực của kết quả FFT, sau đó nó là phổ biến hơn nhiều để trích xuất độ lớn của FFT bằng cách sử dụng abs (F). Điều này tương đương với sqrt (real (F).^2 + imag (F).^2). **

Hầu hết thời gian bạn sẽ muốn sử dụng một NFFT ngắn hơn. Điều này có thể là do bạn có thể chạy phân tích trong một hệ thống thời gian thực, hoặc bởi vì bạn muốn trung bình kết quả của nhiều FFT với nhau để có ý tưởng về phổ trung bình cho tín hiệu thay đổi thời gian, hoặc bởi vì bạn muốn so sánh phổ các tín hiệu có thời lượng khác nhau mà không lãng phí thông tin. Chỉ cần sử dụng lệnh fft với giá trị của NFFT < số lượng các phần tử trong tín hiệu của bạn sẽ dẫn đến một fft được tính toán từ các điểm NFFT cuối cùng của tín hiệu. Đây là một chút lãng phí.

Ví dụ sau có liên quan nhiều hơn đến ứng dụng hữu ích. Nó cho thấy làm thế nào bạn sẽ chia một tín hiệu thành các khối và sau đó xử lý từng khối và trung bình kết quả:

%//These are the user parameters 
durT = 1; 
fs = 22050; 
NFFT = 2048; 
sigFreq = 300; 

%//Calculate time axis 
dt = 1/fs; 
tAxis = dt:dt:(durT-dt); 

%//Calculate frequency axis 
df = fs/NFFT; 
fAxis = 0:df:(fs-df); 

%//Calculate time domain signal 
x = cos( 2*pi*sigFreq*tAxis ); 

%//Buffer it and window 
win = hamming(NFFT);%//chose window type based on your application 
x = buffer(x, NFFT, NFFT/2); %// 50% overlap between frames in this instance 
x = x(:, 2:end-1); %//optional step to remove zero padded frames 
x = ( x' * diag(win) )'; %//efficiently window each frame using matrix algebra 

%// Calculate mean FFT 
F = abs( fft(x, NFFT)/sum(win) ); 
F = mean(F,2); 

subplot(2,1,1); 
plot( fAxis, 2*F ) 
xlim([0 2*sigFreq]) 
title('single sided spectrum') 

subplot(2,1,2); 
plot( fAxis-fs/2, fftshift(F) ) 
xlim([-2*sigFreq 2*sigFreq]) 
title('whole fft-shifted spectrum') 

tôi sử dụng một cửa sổ Hamming trong ví dụ trên. Cửa sổ bạn chọn phải phù hợp với đơn đăng ký http://en.wikipedia.org/wiki/Window_function

Số tiền trùng lặp bạn chọn sẽ phụ thuộc phần nào vào loại cửa sổ bạn sử dụng. Trong ví dụ trên, cửa sổ Hamming cân trọng lượng các mẫu trong mỗi vùng đệm về phía không trung tâm của mỗi khung. Để sử dụng tất cả thông tin trong tín hiệu đầu vào, điều quan trọng là sử dụng một số chồng chéo. Tuy nhiên, nếu bạn chỉ sử dụng một cửa sổ hình chữ nhật đơn giản, sự chồng chéo trở thành vô nghĩa vì tất cả các mẫu đều có trọng số như nhau. Bạn càng sử dụng nhiều chồng chéo thì cần phải xử lý nhiều hơn để tính phổ trung bình.

Hy vọng điều này sẽ giúp bạn hiểu.

Answer 2

Kết quả của bạn hoàn toàn đúng. Tính toán trục tần số của bạn cũng đúng. Vấn đề nằm trên thang y trục. Khi bạn sử dụng hàm xlims, matlab sẽ tự động tính toán lại tỷ lệ y để bạn có thể xem dữ liệu "có ý nghĩa". Khi các đỉnh cosin nằm ngoài giới hạn bạn chọn (khi f> 500Hz), không có đỉnh để hiển thị, do đó thang đo được tính dựa trên một số tạp nhiễu nhỏ (ở đây tại máy tính của tôi, với MATLAB 2011a, thang y là 10 -16)

Thay đổi giới hạn thực sự là phương pháp chính xác, bởi vì nếu bạn không thay đổi nó, bạn không thể thấy các đỉnh trên phổ tần số.

Một điều tôi nhận thấy, tuy nhiên. Có lý do nào để bạn vạch ra phần thực sự của biến đổi không? Thông thường, nó là abs(F) được vẽ và không phải là phần thực.

chỉnh sửa: Trên thực tế, bạn là trục tần số chỉ đúng vì df, trong trường hợp này, là 1. Dòng fax là đúng, nhưng tính toán df thì không.

FFT tính điểm N từ -Fs/2 đến Fs/2. Vì vậy, N điểm trên một loạt các Fs sản lượng một df của Fs/N. Như N/thời gian = Fs => thời gian = N/Fs. Thay thế nó trên biểu thức của df bạn đã sử dụng: your_df = Fs/N * (N/Fs) = (Fs/N)^2. Như Fs/N = 1, kết quả cuối cùng là đúng: P