texts: Duality of the continuous-time Fourier transform - derivation and notation

texts:
Suppose we have the Fourier transform pair $x(t)$ and $X(\omega)$ such that
 $$X(\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \mathrm{d}t$$

The duality property states that $X(t)$ and $2\pi x(-\omega)$ constitute a Fourier transform pair. I was trying to prove this statement when I ran into the following problem: the Fourier transform of a signal is sometimes denoted by $X(j\omega)$, such that $$X(j\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \mathrm{d}t$$

Using this notation, how do we even state the duality property? Note that evaluating $X(t)$ in this case would effectively replace $j\omega$ (as opposed to just $\omega$ in the previous case) by $t$ which is clearly wrong. Would we have to write $X(jt)$ for consistency?

Returning to the question of deriving the duality property and using the notation I used in the beginning, what is the flaw in the following approach?

Since $$X(\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \mathrm{d}t$$
we let $t$ be dummy variable $u$, and evaluate the function $X(\cdot)$ at $t$. We then have $$X(t) = \int_{-\infty}^{\infty} x(u) e^{-j ut} \mathrm{d}u$$
With the change of variables $\omega = -u$ we get $$X(t) = -\int_{-\infty}^{\infty} x(-\omega) e^{j\omega t} \mathrm{d}\omega$$ $$X(t) = -\frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi x(-\omega) e^{j\omega t} \mathrm{d}\omega$$

Which is almost the correct result except for the additional minus sign. Where did I go wrong? Thank you in advance.

I would also like to add that I feel uncomfortable about simply changing $t$ and $\omega$ in the functions' argument because they have different dimensions.


texts: Can we use Tiva C series TM4C123GXL for signal processing experiment purposes?

texts: I've newly bought texas instruments Tiva C series TM4C123GXL processing board for developing DSP application project. Is that processing board is good for that purpose? Is that well work with Matlab? Or C2000 Piccolo TMS320F28075 32-bit microcontroller is good for my purpose ?


texts:
Is it possible to calculate the SNR from finite duration time-domain samples that contain the received signal, $y(t)$, where:
$$    y(t) = x(t) + n(t) $$
and $x(t) =$ transmitted signal, and $n(t) =$ AWGN noise
The signal $y(t)$ is complex baseband, filtered and resampled at Nyquist. The question here is, can the $\mathbb{E}\big[n(t)\big]$ be calculated from the finite number of samples that contain the energy of the signal?
note: this question addresses how to do it if you have periods of time where the signal is not transmitting: SNR calculation


texts: What are the fundamental reasons that suggest us to minimise transmission power for energy-efficient radio link design?

texts:
My question in short:

Under which scenarios minimising transmission power is beneficial in terms of transmitter energy-efficiency? is there any fundamental reason to keep a low-power oriented design strategy?

In particular, small, low-cost transceivers (such as those used for M2M/IoT communications) are designed according to this strategy. I'd like to understand why.

Background

It is quite natural to think that energy efficiency in radio link design is directly linked to a low transmission power. This, however, cannot be concluded from a simple analysis (see baseline model below). For instance, we may feel tempted to decrease the data rate in order to reduce the required transmission power. This in turn increases the transmission time, and the total energy consumed does not change.

Naturally, this question requires a look to the whole transmitter circuitry (which is a little outside my field). After some research, I noted that there must be a trade off regarding clock rate (high-precision and high-rate is payed with extra power) and amplifier efficiency (lower input powers are usually less efficient). But this suggests that minimising power is not always the best option.

The baseline model (optional)

Consider a wireless radio system where a transmitter $A$ wants to communicate with a receiver $B$, at distance $d$. The total attenuation (due to to signal propagation) can be simply modelled as $L=d^\alpha$, where typically $2 \leq \alpha \leq 4$.

The carrier-to-noise density ratio $\frac{C}{N_0}$ at $B$'s receiver input is

\begin{equation}
 \frac{C}{N_0} = \frac{P_t G_A G_B}{d^{\alpha}N_0},
\end{equation}

where $P_t$ is the power at the output of $A$. When designing these systems, there is usually a performance constraint expressed as a maximum tolerated bit error rate, which turns into a minimum $E_b/N_0 = \gamma_0$ requirement.

So in order to achieve this $\gamma_0$, it seems that there are only two system design parameters to play with (assuming all other parameters are fixed): transmission power and data rate, i.e.

\begin{equation}
 \gamma_0 \propto \Big( \frac{P_t}{R_b} \Big).
\end{equation}

Thus, we observe that increasing the transmission power has the same effect that decreasing the data rate,  in terms of $E_b/N_0$.

Now consider that the total amount of data to be delivered is $L_b$. The total energy consumed by $A$ to deliver $L_b$ bits is therefore

\begin{equation}
 E_t \propto  P_t\frac{L}{R_b}=\gamma_0 N_0 L,
\end{equation}

which does not depend on the operating power neither on the data rate.


texts: BFSK modulation what signal the transmitter actually sends?

texts:
I am trying to understand how modulation with carriers works(transmission part). I will express my question in terms of BFSK for simplicity.

Let us condiser BFSK modulation with symbols {1,2}.

Then for a single symbol period:

the transmitted waveform for symbol 1 is

$$s_1(t)=\sqrt{P\over T}e^{j2\pi f_1 t}$$

the transmitted waveform for symbol 2 is

$$s_2(t)=\sqrt{P\over T}e^{j2\pi f_2 t}$$

where $P , T$ denote signal strength and period respectively and $f_1 , f_2$ denote carrier frequencies for symbols $1 , 2$ respectively .

We also know that $s_1(t)$ has an equivalent
$$s_1'(t)=\Re{(s_1(t))}=\Re{(\sqrt{P\over T}e^{j2\pi f_1 t}})=\sqrt{P\over T}cos(2 \pi f_1 t)$$

Same thing for $s_2(t)$

$$s_2'(t)=\Re{(s_2(t))}=\Re{(\sqrt{P\over T}e^{j2\pi f_2 t}})=\sqrt{P\over T}cos(2 \pi f_2 t)$$

Question 1: When I send symbol $1$ which signal is actually transmitted from the transmitter, $s_1(t)$ or $s_1'(t)$?

Question 2: For signals $s_1(t)$ and $s_1'(t)$ which one of two is the low-pass equivalent and which one is the highpass equivalent?


texts:
I was wondering by what logic or mathematics the frequency-ration $$ \frac{f_{lo}}{f_s} $$ in the sine-functions comes into existence?

Source


texts:
Is there any difference between the two?

When we say discrete signal, does it only imply that the x-axis or time is discrete? Or it can also mean discrete amplitudes on the y-axis?

Or is the signal is called digital only when both of its axes are discrete?


texts: Estimating Real-Time Vehicle Velocity from low-frequency GPS, accelerometer, gyro and magnetometer in Smartphone

texts:
I am trying to write a system that would be able to estimate the velocity of a traveling vehicle (usually a car) using a smartphone. Most smartphones have the next sensors:

1Hz Low-Frequency GPS
10-100Hz Accelerometer
10-100Hz Gyroscope
10-100Hz Magnetometer

I would also impose a limitation, that device must be in either portrait or landscape position, relative to vehicle's direction of travel.

Now, the GPS will tell me speed and position once per second, but I would like to estimate the velocity in time between GPS data comes in, up to 10-100 times per second, for which I would like to use other sensors. The speed should be accurately represented, even if the car is decelerating.

For now, I've implemented Kalman filter which fuses together GPS, accelerometer and gyroscope and estimates the velocity. But the biggest problem becomes the deceleration, as in all usual references I could find, only the absolute velocity from all axes is taken into consideration.

How to approach this problem? Am I even correct in using Kalman? How to get correct deceleration? What if vehicle is cornering?

Thank you for your help!


texts: How to estimate a transfer function from a magnitude-only frequency response?

texts:
Given an arbitrary frequency response, what signal processing methods might exist that could guess, estimate or determine a transfer function (pole and zero constellation) which gives a "reasonably good" approximation (for some given estimation quality criteria) to that given frequency response?  What means exist to estimate the number of poles and zeros required for a given transfer function plus a given approximation error allowance?  Or how could one determine these constraints can't be met, if possible?

If the given frequency response was actually produced by a known transfer function, will any of these methods converge on that original transfer function?  How about if the given frequency response were subject to (assumed Gaussian) measurement errors?

Assume working in the Z-plane with sampled spectrum, although continuous domain answers might also be interesting.

Added: Are the solution methods any different if only the magnitude of the frequency response is given (e.g. a solution with any phase response is allowed)?

Added: The latter problem is what I'm most interested in, given a known magnitude response around the unit circle, but unknown/unmeasured phase response, can the measured system be estimated, and if so under what conditions?


texts: Magnitude spectrum for discrete-time signals

texts:
I am doing the following question (it is not homework, I am preparing for an exam).

  The time discrete signals $x_1(n)$ and $x_2(n)$ is created by sampling the continuous signal $x_a(t) = \cos (2 \pi 300t) + \cos (2 \pi 600 t)$ with sampling frequency $F_s = 1000\textrm{ Hz}$, $x_1$ with anti aliasing filter and $x_2$ without aliasing filter.

  a) Calculate $X_a(f)$ for the time continuous signal $x_a(t)$.
  b) Calculate the magnitude spectrum for the time discrete signals $x_1(n)$ and $x_2(n)$ i.e $\lvert X_1(f)\rvert$ and $\lvert X_2(f)\rvert$, with $-1<f<1$.

I have done a), I used the table for Fourier Transform to conclude that $$X_a(f) = \frac{1}{2}\Big(\delta(f+300) + \delta(f-300) + \delta(f+600) + \delta(f-600)\Big).$$

But I can not solve b), I believe that I can use the information/answer from a) that is why I provided it.


texts: Lowpass Butterworth filter design (own poles' calculations)

texts:
I'm trying to calculate coefficients for my Butterworth filter of my own. I've found this article that looks really reliable. But I'm trying to implement it in MATLAB and I'm getting bandstop filter.
Please find attached matlab script below:

clc, clear all, close all;
%DIGITAL BUTTERWORTH FILTER DESING
%Butterworth LP with -3dB at fc and N order and G gain
fs = 10; %frequency sampling
fc = 1;   %cut off frequency
N = 2;      %filter order
G = 1;      %filter gain
wc = 2*pi*fc;

%analog filter design
%poles and zeros
k = [1:1:N];
%poles searching
%sk = wc * exp(((j*pi)/(2*N))*(2.*k+N-1))
sk = wc * exp(j*(pi/2)) * exp(j*(2.*k+1)*pi/2*N);
z = 0;
G = 1;
[NUM,DEN] = zp2tf(z,sk,G);
%no zeros for LP
Hs = tf(1,DEN)

%filter analysis
figure(1);
subplot(3,1,1)
pzmap(Hs);
subplot(3,1,2)
impulseplot(Hs);
subplot(3,1,3)
stepplot(Hs);
fvtool(1,DEN);

I've calculated poles according to equation $10.3.47$ in the article and created transfer function with this poles as a parameters. But fvtool shows bandstop filter instead of lowpass (as described in article). In pzmap I'm getting also pole on right side of $s$-plane so it seems to be unstable. I dont really know how to interpretet this results.

Maks


texts: How does an Image filter really sweep an Image?

texts:
Suppose, we have the following image:

Now, we want to apply a 5x5 filter to the image.

How would the filter sweep the image?

Like this?

Or, like this?

Or, anything else?


texts:
I am working on getting room impulse response(RIR) and the first order reflections of a room.
For now I am stuck at getting proper RIR.
I followed Joseph's work and implemented on my own.

My setup uses:

microphone - mic in webcam (Logitech C920 HD Pro)
speaker - usb powered speaker connected to 3.5mm laptop audio port (havit brand with squarish shape)
room - almost empty room. dimension about 6.3 x 4.8 m or 250" x 190"
audacity - edit audio obtained
matlab - plotting of the graphs and doing other stuffs

The image below contains four graphs;

First is the recorded sound by Joseph
Second is the impulse response of a room used by Joseph
Third is the recorded sound by me
Fourth is the impulse response I've gotten.

The impulse in second graph is quite visible as compared to the fourth graph.
Here is my question.

Why is the impulse in fourth graph (the one I did) not very visible?
What causes this to happen?

I suspect it is caused by equipment and room as that is the only difference i can think of.

I wish to listen what you guys would say about this. Feedbacks, guidance and improvements will be very much appreciated.

Thanks


texts:
I am working on the project where my goal is to create SR image from 25 LR images that are precisely shifted by 0.2px (both in horizontal and vertical direction ). I also know that the captured object is a back and white grid ( something like chessboard)
I am looking for some precise SR method that will take advantage of the precise shift of the images. The emphasis is on precision not on the speed of the method. I have read a lot of the papers on that topic and found out that there is a lot of methods and it is difficult for me to decide which method suits my task.

If you have experience with SR please suggest some methods that find good for this task.


texts: How to plot Full spectrum with negative frequencies

texts:
I'm trying to do frequency analysis on signal of a rotating shaft. The signal is generated in two measurement planes i,e X-position and Y-Position. I was able to use fft and transform into frequency domain for individual planes. The plot I get for each place is only Half spectra. What I'm looking for is a full spectra which could combine both the planes. I did some research and found the following description online.

Full spectrum: The full spectrum is an additional diagnostic tool and is also called the spectrum of an orbit. It shows the same information as an orbit but in a different format. It helps to determine the degree of ellipticity (or flattening) associated with the various machinery conditions along with the precessional direction for all the frequency components present. To obtain the full spectrum, the orthogonal X and Y transducer signals are fed into the direct and quadrature parts of the FFT input. The positive and negative vibration components for each frequency are obtained. Positive is defined to be the forward precession and the negative component as the reverse precession. These components yield the following ellipticity and precessional information for a given orbit of any particular frequency (1× or 2× or …):

• The sum of two components, forward and reverse, is the length of the orbit major axis. • The difference between the two components is the length of the orbit minor axis. • The larger of the two components, positive or negative, determines the direction of precession that is forward or reverse.

One of the possible applications of full spectrum is analysis of the rotor runout caused by mechanical, electrical or magnetic irregularities. Depending on the periodicity of such irregularities observed by the X–Y proximity probes, different combinations of forward and reverse components are observed. The method forms the basis for many useful machinery diagnostics.

The full spectrum (just like the normal FFT) can be obtained in a steady-state analysis (a single FFT or waterfall) and even in transient analysis, which would then be called the full spectrum cascade (Figure attached).

Could anyone help how the 2 half spectra is converted into full spectra plot. Does anyone have any idea about this? Any hints or suggestions please!!


texts: Digital filter not removing noise at specific frequency Matlab

texts:
Doing some work at the minute on digital filters in matlab, I have a file with artifical noise added (sine wave added at specific frequency). The goal is to filter the signal and get it as close as possible to the clean signal provided.

I've done an FFT and plotted the results and found a very large spike at 29.3Hz which is not present in the clean signal.

I've tried using a notch filter, which I thought would work since it operates at such a specific frequency, however it just seems to attenuate the signal and remove some power but not block it completely. I then added a bandstop filter to try and block any signals in that region and it simply attenuated the signal also. Does anyone have any thoughts? I just seem to be lowering power of the entire signal and not actually removing anything, getting the basic shape of the clean signal but still a lot of noise present after both filters. Thanks!

[b1,a1] = iirnotch((29.3*(2/fs)),0.99999);
IIR1 = filter(b1,a1,ecg58_DC_removed);

FFT_resultFilter1 = (1/length(t))*fft(IIR1);
f=(0:1024)/1024*(200/2);
figure(4)
stem(f, 2*abs(FFT_resultFilter1(1:1025)));
xlabel ('Frequency (Hz)');
ylabel ('Spectral Magnitude');
title('First filter')
grid on

[a2,b2] = butter(2,[29.2 29.4]*2/fs, 'stop');
IIR2 = filter(a2,b2,IIR1);

FFT_resultFilter2 = (1/length(t))*fft(IIR2);
f=(0:1024)/1024*(200/2);
figure(5)
stem(f, 2*abs(FFT_resultFilter2(1:1025)));
xlabel ('Frequency (Hz)');
ylabel ('Spectral Magnitude');
title('First filter')
grid on

figure (6)
plot(t(1:1000), IIR2(1:1000));
xlabel('time (s)')
ylabel('amplitude')
title('two filters');

b1 =

1.0e-04 *

0.1571   -0.1902    0.1571

a1 =

1.0000   -0.0000   -1.0000


texts:
I am trying to understand how to get time as I already have sample rate. I want to draw a graph of sample rate vs time. I have information such as carrier frequency, channel bandwidth. Can that be used to get time.

So far, I have understood this: Samples per symbol and number of symbols for QAM

In order to get the time, I would need to do this:

T = 1/Samples per symbol

Would this be correct.


texts: Is Karplus-Strong synthesis a reliable way to test pitch detection algorithms?

texts:
I'm using Arduino DUE (Cortex-M3 based board).

I want to make a guitar sound using DAC (digital to analog converter).

For instance, the note A2's fundamental frequency is 110 Hz (using A4 440 Hz as the standard).
As you know, guitar generates harmonics like 220, 330, 440 Hz, and more.
More complex, the amplitudes of each sine waves change as time goes by.

This is described with detail at this paper (Figure 1, Electric & Acoustic Guitar Strings: A Recording of Harmonic Content by Ryan Lee)

I'm trying to understand the Karplus–Strong string synthesis.

(I used Matlab and some examples from Mathworks to understand it.)

However, I'm not sure whether using this algorithm will generate a guitar sound like the acoustic guitar.

What I mean is that, I'm trying to test my pitch detection algorithm with a guitar. But taking a guitar all the time for testing is quite bothersome.
So I tried to generate a wave using DAC for faster testing.

When using DAC, will the generated waveform which uses Karplus–Strong algorithm be a reliable input for testing a instrument pitch detection algorithm?

Or is there something else that can be used?


texts: Filtering kinda "half-toned" image for OCR processing

texts:
I have a scanned PDF material to which I want to add hidden text layer, so I could index the document. I used ghostscript black and white tiff output device (tiffg4) to extract pages as tiff images, and here is example of what they look like:

Processing this image with tesseract, does not give good results.
Changing ghostscript output DPI (600, 300, 150, 96) shows that image at 96 DPI gives best result from tesseract but it's still not satisfactory.

Now I thought to ask for advice which filter would enhance this image for OCR processing.

I could use imagemagick, or numpy/scipy/ndimage


texts: DTMF Tone Generation in MATLAB not sounding right

texts:
I'm starting out with some DTMF Tone Generation using Matlab and I wanted to find some theory behind how to generate the various tones for a 4x3 keypad consisting of 0-9, * and #. So far I have the following function:

function y = OutputTone(c)

t = 0:0.1:500;
num = double(c);
%fs = 8000;

if (num == 49) || (num == 50) || (num == 51) % 1, 2, 3
    f1 = 697;
    if (num == 49)
        f2 = 1209;
    elseif (num == 50)
        f2 = 1336;
    elseif (num == 51)
        f2 = 1477;
    end
elseif (num == 52) || (num == 53) || (num == 54) % 4, 5, 6
    f1 = 770;
    if (num == 52)
        f2 = 1209;
    elseif (num == 53)
        f2 = 1336;
    elseif (num == 54)
        f2 = 1477;
    end
elseif (num == 55) || (num == 56) || (num == 57) % 7, 8, 9
    f1 = 852;
    if (num == 55)
        f2 = 1209;
    elseif (num == 56)
        f2 = 1336;
    elseif (num == 57)
        f2 = 1477;
    end
elseif (num == 42) || (num == 48) || (num == 35) % *, 0, #
    f1 = 941;
    if (num == 42)
        f2 = 1209;
    elseif (num == 48)
        f2 = 1336;
    elseif (num == 35)
        f2 = 1477;
    end
end

y = sin(2*pi*f1*t) + sin(2*pi*f2*t); % DTMF Signal in Time Domain

end

As the code describes above, I take in an input character (0-9,*,#), convert it to ASCII and then select the appropriate DTMF frequency pertaining to the key that is pressed. For example, if 1 is input to the function, then the frequencies f1 = 697Hz and f2 = 1209Hz will be selected and applied to the output sinusoid y.

However, when I play back the sound of y, it doesn't sound anything like the key tones I hear on my smartphone or in this link. What is the correct way to represent y?


texts: What does following statement mean: " Simulation was done at the rate of 12 samples per cycle of 50Hz a.c."

texts:
I am very new to the signal processing world and i am currently doing a project of studying frequency components in electrical machines. I am simulating the system in PSCAD. But i need to simulate the model at the rate of 12 samples per cycle of 50 Hz. What does that statement mean?
TIA


texts: How to compute transfer function $G(s) = \exp \left( - \sqrt{s} \right)$ in Matlab / Simulink?

texts:
How to compute transfer function
$$G(s) = \exp \left( - \sqrt{s} \right)$$
in Matlab / Simulink?
I am trying to calculate a PID controller for this function. This function describes heat transfer via conduction. I want so see first the step time response and impulse response for this system. I usually use for other functions the tf function from Matlab, including IODELAY, but in this case i can't compute $\sqrt{s})$.The error is:Undefined function 'sqrt' for input arguments of type 'tf'. Is there any other way to compute this transfer function?

Hi again, i just edited it because i found more information.In my case i guest that T=1.

For Figure 2.21 i have this:


texts: In the literature, I've found that DFS and DFT are one and the same. If they are one and the same why to use two different names for them? If there is really a difference what is it and what is the significance of discrete Fourier series?


texts: How to prove relationship between mutual information and differential entropy?

texts:
I need to prove that mutual information given by

$$I(X;Y)=\int_{x,y}f(x,y) \log_2 \left( \frac{\left (f(x,y)\right)}{f(x) f(y)}\right) \, dx \, dy$$

is equivalent to $I(X;Y)=H(Y) - H(Y|X)$ I am proceeding like this

$$I(X;Y)=\int_{x,y} f(x) \frac{ f\left(x,y \right)}{f(x)} \log_2 (f(Y|X)) \, dx \, dy - \int_{x,y}f(x,y) \log_2(f(y)) \, dx \, dy$$

where first term gives me $H(X|Y)$ and second term gives me

$$\int_{x,y}f(y) f(x|y) \log_2(f(y)) \, dx \, dy $$

and I am taking

$$\int_{x} f\left (x|y \right) dx =1$$

I mean to say if this equation is right what does this physically symbolize that given $y$ there must be some $x$ that means given output there must be some input with probability $1$? Am I doing the derivation right ?


texts: Does constant Q transform have linearity property in the transformed domain?

texts:
With FFT, I sometimes take advantage of the fact that I can pre-calculate a signals Fourier Transform, and then you can add the noise in the spectral domain:

$$ \mathscr{F} (x[n] + h[n]) = \mathscr{F}(x[n]) + \mathscr{F}(h[n]) $$

This is useful when we want to discard of the original signal, but we still want to noise it with a spectrum. If we know the spectrum of the noise beforehand, this is also faster computationally (just adding together the two spectra).

Does the same relation hold for constant Q transforms?


texts: Question about relationship between aliasing and matched filter

texts:
while I was reading J. Barry's textbook (https://www.springer.com/gp/book/9780792375487), page 308 says the following:
"... symbol-rate sampling at the matched filter output is generally at less than the Nyquist rate, so aliasing of both noise and signal is inherent in this sampling. This aliasing will not compromise the performance of the receiver as long as the filter before the sampling is a matched filter."
Here, the context of this text is about how whitened matched filter (WMF) can provide sufficient statistics for ML detector. But I do not understand why he mentions that aliasing of signal and noise does NOT degrade the system performance as long as the matched filter is deployed.
My suspicion is that noise is whitened after WMF, and hence the aliasing will simply increase the variance of the noise but does not change the autocorrelation between noise samples (here, since the noise is whitened, its autocorrelation should be 0). But the following thought was that the ISI generated by the aliasing will decrease the system performance at the end...
How am I supposed to understand this text?


texts: Generate time-domain random signal from PSD

texts:
Given an analytical description of the PSD, for example (MATLAB "pseudocode"):

N = 20400;
bw = 2 * 2040 / N; % double-sided spectrum, 20% occupied
PSD = zeros(N,1);
P0 = sigma * sqrt(N/bw); % sigma^2 = (P0^2/N) * bw
PSD(1:2040) = P0;
PSD(18362:end) = P0;

what is the "correct" way to generate a time-domain signal $x[n]$ which is a realization of this PSD, i.e. which has variance of $\sigma^2$ and a DFT where all bins $2041, \dots, 18361$ are zero?

I found this: http://www.mathworks.com/matlabcentral/newsreader/view_thread/264846. However, I am confused since there only the phase is taken randomly but not the magnitude which would result in a perfect "brick-wall" spectrum when taking the DFT of each realization (because the magnitude and magnitude-squared is not random). However, this is only the case if the number of samples approaches infinity.

Important: I do not want to generate iid. white Gaussian samples and filter them - I would like to generate the random signal in DFT domain and then use the inverse DFT to generate the time domain signal.


texts: Discard circular/blob like structures from image

texts:
I have a binary image that looks like this:

And I would like to get rid off all the small structures and just keep what's on red:

Image erosion is an option that works well to remove the very small groups of pixels without affecting the red regions. Is there any other method that can be applied to keep only the elongated structures on my image?


texts: issue in DFT of OFDM by blocks and N-by-N

texts:
I'm trying to study and understand the OFDM, and during coding the OFDM in MATLAB, I faced this issue:

suppose we have the Numsymbol = 256, and dataqpsk our used data, The issue which I face that when I perform iFFT block by block, that seems to be ok, but when trying to make it N-by-N .. that gives a problem.

here is the steps. now I will make the iFFT by blocks :

clear all; clc;

N_carrier = 256;  % number of carrier
N_symbol = 256;  %number of symbols
N_FFT = 256;  %length of FFT
N = N_carrier;
L_cp = N/4;   % the CP

dataqpsk = rand(N_carrier,N_symbol);  %Example for modulated data
size(dataqpsk);
for isij=1: N_symbol
    % input modulated data
    dataqpsk1=dataqpsk(:,isij);
    size(dataqpsk1)
    D_map = dataqpsk1;

    %% doing such mapping
    D_IFFT_fore=zeros(N_FFT,1);
    D_IFFT_fore(1:(N/2)) = [D_map(((N/2)+1):N)];
    D_IFFT_fore((N_FFT-((N/2)-1)):N_FFT) = [D_map(1:(N/2))];
    D_IFFT_aft=N_FFT.*ifft(D_IFFT_fore,N_FFT);%×÷ifft

    %add GI and interpolation.
    D_cp=[D_IFFT_aft(N_FFT-L_cp+1:N_FFT);D_IFFT_aft];
end

But now when I try to make it as DFT N-by_N I don't get, as below:

clear all; clc;

N_carrier = 256;  % number of carrier
N_symbol = 256;  %number of symbols
N_FFT = 256;
N = N_carrier;
L_cp = N/4;  %CP

dataqpsk = rand(N_carrier,N_symbol);
size(dataqpsk);
D_map = dataqpsk;
D_IFFT_fore=zeros(N_FFT,N_FFT);
D_IFFT_fore(1:(N/2)) = [D_map(((N/2)+1):N)];
D_IFFT_fore((N_FFT-((N/2)-1)):N_FFT) = [D_map(1:(N/2))];
D_IFFT_aft = N_FFT.*ifft(D_IFFT_fore,N_FFT); % ifft
D_cp(1:L_cp,1:N_FFT)=D_IFFT_aft(N_FFT-L_cp+1:N_FFT,1:N_FFT);

The problem is when trying get the results of DFT N-by_N, I get from the step of D_IFFT_fore(1:(N/2)) = [D_map(((N/2)+1):N)];  a matrix of Zeros !!

please, any help?

thank you

Thank you Dear msm, and here is the complete code:

clear all; clc;

N_carrier = 256;  % number of carrier
N_symbol = 256;  %number of symbols
N_FFT = 256;
N = N_carrier;
L_cp = N/4;  %CP

dataqpsk = rand(N_carrier,N_symbol);
size(dataqpsk);
D_map = dataqpsk;
D_IFFT_fore=zeros(N_FFT,N_FFT);
D_IFFT_fore(1:(N/2),:) = [D_map(((N/2)+1):N,:)];
D_IFFT_fore((N_FFT-((N/2)-1)):N_FFT,:) = [D_map(1:(N/2),:)];
D_IFFT_aft = N_FFT.*ifft(D_IFFT_fore,N_FFT); % ifft
D_cp(1:L_cp,1:N_FFT)=D_IFFT_aft(N_FFT-L_cp+1:N_FFT,1:N_FFT);

thanks again


texts: Applying a high-pass filter on constant signal

texts:
I recently implemented a Butterworth high-pass filter (2nd order). Everything seems to work fine, except that I have a question when applying such a filter to a constant signal. As there are no high-frequency components in a constant signal, I would expect the filter to yield a constant 0 signal.

In the plots below are my results. There seems to be ``ripple'' in the first couple of frames before the high-pass filtered signal (correctly) converges to 0.

Is this a logical result from applying a Butterworth high-pass filter on a constant signal, or might there be a bug in my code?

EDIT: I've created the same filter in Octave, resulting in the same output. This indicates that we're indeed looking at the step response, as Paul R. indicates.

Octave output


texts: mne raw plot looks like vertical lines instead of the usual signal

texts:
I have a data from a single eeg channel. When I plot it with plt, it looks like a regular signal.

However, when I'm trying to plot with mne raw object's method plot, it looks very strange, like vertical lines:

Does anyone know why it looks like that?
Thanks!


texts: All-pass filter usage in automotive room correction

texts:
This may not be the right place to ask this question but I'll give it a shot anyways. Maybe someone here is an audiophile.

I'm using a DSP called Audison Bit One HD which uses FIR filters for crossover and EQ. They recently included all-pass filters which allows me to adjust the phase of the crossover point. I'm wondering if anyone here knows how this will benefit for room correction.

I know what an all-pass filter does but I'm not sure how it could be useful with FIR filters which is a linear phase filter. Why would I want to adjust the phase of the crossover point for?

The manufacturer is no help.


texts: How to estimate the poles of a linear combination of exponentials?

texts:
Consider a signal that can be expressed as a linear combinations of exponentials:

$x(n) = \sum_{k=1}^K c_k z_k^n$

I would like to estimate the poles $z_k$ and I was told that I could estimate them by calculating the SVD of the Hankel matrix of the observations $\mathbf{H}$:

$\mathbf{H} = \left(\begin{array}{cccc}x(1) & x(2) & x(3) & \dots \\
x(2) & x(3) & x(4) & \dots \\
x(3) & x(4) & x(5) & \dots \\
\vdots & \vdots & \vdots &  \\\end{array}\right)$

But how exactly do I utilize the SVD? The poles will be the singular values? And is it possible to find $K$?

Thanks in advance!


texts:
I am trying to read OFDM synchronization paper by Schmidl and Cox. The math is slightly overwhelming.

What is the formula that they use in implementation that distinguishes between coarse estimation and fine estimation?

Secondly once fine frequency offset is done, how do you determine that you need to do coarse frequency synchronization or not?  What is it that you will check to determine it? If anybody could explain with a set of equations that would  be very helpful.


texts: Welch's overlapped method for amplitude by frequency estimation

texts:
pwelch  is a useful matlab function for computing real/magnitudes for power using WOSA method, but I'd like to apply the same thinking in order to output complex pairs for amplitude (not power) spectra.

Does anyone know if this is possible? Or if there is an existing script that I might be able to modify?

Thanks!


texts: Resampling with factor linearly moving from one value to another: calculate expected target length

texts:
I am resampling an audio (1D) signal, using a resampling factor that moves linearly from startFactor to endFactor; say the input signal has a length of 4194304 sample frames, and startFactor is 1.0 and endFactor is 4.0, for example. Since the DSP blocks reads the resampling factors from a signal (here a linear ramp) at one sample per resampled output sample, I need to know the target length in advance to generate a linear movement.
My first intuition was that this is the geometric mean of startFactor and endFactor, but that is unfortunately not the case. For the given example values, that would be 2.0, and my ramp would move from startFactor to endFactor over 8388608 frames. However, the resulting output "overshoots" and the target is 9656320 sample frames instead (+15%). Likewise, if I change endFactor to 0.25, the geometric idea would estimate the target to have length 2097152, but the actual length is 2176638 (+4% overshoot).
Those discrepancies tell me it's not just floating point noise, but the formula has to be different.
I know I can calculate the target length by integrating the factor signal, but I need the target length to create the factor signal, so it's a catch-22.

I can also iterate to "find" the value, but I need an analytical solution. For example, the correct length factor for 1.0 and 4.0 is ca. 2.16356, so that both the ramp and target signal have length 9074628. For 1.0 and 2.0 the length factor is ca. 1.44236, for 1.0 and 0.5, the length factor is ca. 0.721344, and for 1.0 and 0.25, the length factor is ca. 0.541011.


texts:
In this paper, the following prewhitening filter is described:

$$
C(z) = \sum_{k=0}^n c_{k}z^{-k}
$$

where $n$ and $c_k$ are known.

The paper also describes the values $C(\lambda_{k})$, with $\lambda_{k}$ complex number, as the above filter "evaluated at the mode $\lambda_{k}$".

Does anybody understand what the $C(\lambda_{k})$ values represent (what do they mean by a filter mode)? Are they just $C(z)$ computed by replacing $z$ with $\lambda_{k}$?


texts: Variations in MATLAB plot each time the program is executed

texts:
Every time I run the program my PLOT in MATLAB gets changed to some extent (some times to larger extent too).
I understood that the channel is random and this might happen due to random distributions used to model the channel. Is there any way to solve this issue, so that every time I get the same values and plot.


texts: Two signals with the same frequency in the same area can interfere?

texts: If one signal is transmitting form transmitter $TX_1$ with some frequency $X$ and some other transmitter $TX_2$ is sending some other signal with same frequency $X$ will there be any clash or interference between these two signals since these two transmitters $TX_1$ and $TX_2$ are located in the same area?


texts: With the same channel S/N, the normalized mean square error differs

texts:
I am working on a research work on transmitting ECG signals over wireless body area network. The signal is affected by noise and small scale fading. The normalized mean square error(NMSE) is used to estimate the quality of the reconstructed signal at the receiver. When I run the same m file in matlab several times, I obtain different values for the NMSE, although the same channel signal to noise ratio is used.

I think this is reasonable, do you agree?
If so, how can I obtain a single value to represent the quality of the reconstructed signal at the receiver?
I want to plot a curve with the S/N at the $x$-axis and NMSE at the $y$-axis, how can I plot this curve if the value of NMSE differ with different runs of the same m file with the same value of S/N?


texts: Calculating output of a system given impulse response and input

texts:
I'm trying to calculate the output given an impulse response $h(t)=e^{-t}u(t)$ and input $x(t)=\cos(2\pi t)$.

I know I need to use the convolution but given that cosine is periodic, I don't see how I would get an answer to converge. Am I thinking about this the right way or missing something?


texts: How to design a delay and sum beam former with fractional delay less than one sample?

texts: I am supposed to design a delay and sum beam former for a short audio. But the sampling frequency is only 8kHz and after times the delay, the fractional delay is less then one sample, and after the "floor" function in Matlab, all the fractional delay integer will be 0 and the results are neither fractional delayed nor time aligned. Exponential weights cannot be used in this case for the reason no imaginary part should be in the filtered signal.


texts: I know that we can design a linear phase FIR filter by choosing filter coefficients symmetric or antisymmetric, what is the intuitive idea behind that?


texts:
If the 16-point DFTs of $x[n]$ and $y[n]$ are given as
$Y[k]=\begin{cases}2X[k], & k=0,2,4,...,14 \\ 0, & k=1,3,5,...,15\end{cases}$, where $x[n],y[n]=0, \forall n<0, n>15$, how can I write $y[n]$ in terms of $x[n]$?

If I take $H[k]=1+(-1)^k$ and multiply it with $X[k]$, I should get $Y[k]$. But when I take its IDFT, it turns out to be zero. How to solve this?


texts: I am reading about different integral transforms. Fourier transform being a subject I will soon take a class on. Looking at the wavelet transform there is one term that I do not understand, that is the noiselet. What would a noiselet look like and how is it defined? Is it some basis of the wavelet transform? It is hard to find information about this subject. Thanks!


texts: Call recognition using Digital Signal Processing

texts:
We are looking for writing applications which would recognize Bird calls and Frog calls using Sound processing.

EDIT: I want to identify the species of bird / frog and do not want to differentiate whether the call is from a bird or a frog. There are about 1300+ bird species in India and about 300+ frog species which we are looking to identify based on the call.

I have a background of iOS application development and C / C++ programming but have never worked on Digital Signal Processing.

We do have few researchers who would collect multiple calls of birds and frogs and categorize it as per species. These calls can be used as reference purpose (?) to identify the bird or frog species by recording the call and comparing this call in user's device.

I do not know the technical feasibility of it, but I would like to know:

Are there any open source solutions which can be extended / used to achieve this?
Is this really possible to achieve with some 80% to 90% accuracy?
If open source solution does not exist, how to go ahead with the implementation?

I understand that the scope of this question is too broad and wide which should be ideally avoided in StackExchange. But would like to have your opinion.


texts: understanding image error metrics, mean squared error

texts:
i am trying to understand the error metrics for images. i had no statistics class or anything and i'm having a bit of problem with even the simplest things.

i have started reading this http://unitstep.net/downloads/school/image-qa-presentation/

just for starters, it says that mean squared error is (page 3)

sum_over_i_j( |ImageA(i,j)-ImageB(i,j)|^2 ) / amount_of_pixels

i understand that this is trying to be something like euclidean distance between two points, but why is there no square root then?

then when you go to the next page of the presentation it tells you that this is a special case of minkowski distance and there IS the square root included in case of mean squared error!

and then, shouldn't we calculate distance of two points by summing the powers of the differences of their dimension components? (`[ (x1-x2)^2 + (y1-y2)^2 + ...] ) and only then square root the sum? then do it for all the pixels and average it.

while on the page 4 it says that we power, sum and square root the pixels first and then do this for all the channels (dimensions) separately?

i don't understand :(


texts: Closed form expression for $\sum_{k=0}^{n} \alpha^{-k}$?

texts: I have seen this sum a couple of times now in the textbook and in the problems book but I don't know how to solve it. I know this might be more suited for math.stackexchange but since this is a fairly common sum and is related I'm pretty sure someone can help. Thank you!


texts:
I have PPG signal that should be normalized, with offline processing I'd just do $z$-normalization,
$$
\displaystyle \frac{x- \mu_x}{\sigma_x} \quad \text{where } x \text{ is the signal.}
$$
The problem is that I get the data in chunks and I'm not sure what would be the effect on the signal when normalizing only chunks of the data at any given point. Any ideas or suggestions on how to approach this would be highly appreciated.

Thanks,
Jack


texts: CIC SINC3 Filter with unsigned Integrators

texts:
I want to use the AD7401 Sigma delta modulator connected to my FPGA.

The Datasheet pg.15 has some demo verilog code which I understand alright except one thing:

They use unsigned integers for all of the registers. Which means especially for the Integrators(that only get 0 or 1 added) that they will definately overflow at some point.
The Integrators of CIC Filters are designed to overflow, but all the sources I found explicitly say that this requires twos complement to work properly.
Am I missing something, or why does the Decimator work with unsigned Integrators aswell(that dont wrap around "nicely").

Also am I right that switching all the registers to signed, and mapping the IO Input (0,1) -> to First Integrator(-1,1) will fix this and result in a signed output.


texts:
I have a C# desktop application.

I am using a wrapper to OpenCV which is EMGU.

I am just performing a basic motion detection operation.

I am using component labelling and 2 frames differencing to determine motion changes.  I am also using an averaging background approach to update the background which i use to 'subtract' from the current and previous frames.

All this works well.

Now consider this image:

Because this bush is close to the camera, movement will be detected by my current methods when the wind blows.

Also, consider this image:

Here we see a cobweb to the right of the image. Again when the wind blows I get motion detection.

My possible solutions:

Look for shapes. In the case of elongated shapes like the cobweb where the height is greater than the width (to a certain scale) then ignore that movement.
But, I found I am missing motion when someone walks on the pavement at the back of the image from across the road.

Erosion.
Again I am missing genuine motion like people and bikes.

Compare Histograms
I have found it is not reliable.  A histogram will only measure the number of different colours (in its own bin).  It is possible for the colours to change due to the saturation and lighting or/and the same numbers of colours are found because another object has passed by.

So, I thought about using atrophy. I did post an earlier question to this here:

Entropy of an image

but as the person who so kindly explained entropy suggested that I should post as a specific/separate question here...

I want to explore whether the frequency or/and entropy would help me here?  Should I look at spatial or frequency distribution?  Would frequency be best to look into if I want to avoid colour changes?

Would either approach be faster than the other as I want it real-time as possible.

Any eduction for me here would be great.


texts:
I'm driving a white noise signal $e(n)$ through an AR random process:

$$ x(n)=e(n)+\sum\limits_{k=1}^N \alpha_k x(n-k) $$

I got five values of $\alpha_k$ and plotted its power spectrum in Matlab for $\omega$ from $0$ to $\pi$.

Then I generated the realizations of $x(n)$ with 10000 samples

Now I'm trying to get the estimate of the spectrum from the expectation data using the periodogram formula

$$ S_{per}(\omega)= \frac{1}{K} \left|\sum\limits_{k=0}^{K-1}x[k]e^{-j\omega k} \right|^2 $$

with $K=256$ the number of first values of $x(n)$. When I run my code what I get is a vector with the same value in every position, so I'm plotting a straight line. I think that just applying the formula I should get the estimate, so what am I doing wrong? :-S

This is the section of the code I'm using to find the estimates

var_e=0.01;                                  % noise variance
b_1=1.4261;                                  % coefficients
b_2=-.7634;
b_3=-.9002;
b_4=1.2548;
b_5=-.5707;
K=256;
wn=rand(1,10000)-0.5;                        % white noise
samp=1:1:1e4;
T=10000;
x=zeros(1,length(samp));
% here I generate the expectations
for n=6:T+5
    x(n)=b_1*x(n-1)+b_2*x(n-2)+b_3*x(n-3)+b_4*x(n-4)+b_5*x(n-5)+wn(n-5);
end
x=x(6:end);
figure
plot(samp,x)
    for i=1:length(w)                         % w is the frequency, has 1000
        for k=1:K
        per2(i)=sum(x(k)*exp(-j*k*i));        % sum of exponentials
        pe2(i)=per2(i)*conj(per2(i));         % squared
        end
    end

Thanks in advance for your help


texts: how to correct minBLEP sawtooth oscillator DC bias

texts:
In my implementation of minBLEP, overlapping minBLEPs result in a DC bias as the frequency increases. This is because each minBLEP correction (i.e. the minBLEP minus the aliased step) has a DC bias.

Here's a frequency sweep from 100Hz to 15kHz. The second signal shows the number of active BLEPs.

Interestingly, the minBLEP paper (http://www.cs.cmu.edu/~eli/papers/icmc01-hardsync.pdf) acts as if the DC woes are all solved by switching from integrating impulses to placing BLEPs. There's no mention of correcting DC when lots of BLEPs overlap.

Should I just add a DC-blocker at the end? Is that what is generally done?

Also, why is it important that the BLEPs are minimum phase? Using linear phase BLEPs would alleviate the DC bias, at the cost of delaying the oscillator output by half the length of the BLEP.


texts: Why is the Fourier transform so important?

texts:
Everyone discusses the Fourier transform when discussing signal processing. Why is it so important to signal processing and what does it tell us about the signal?

Does it only apply to digital signal processing or does it apply to analog signals as well?


texts: Band-edge filtering for an RRC pulse-shaped manchester/biphase-encoded bitstream

texts:
Band-edge (derivative matched) filters can be used to generate an error signal for frequency locked loops (real part) or symbol timing loops (imaginary part). The derivation of the band-edge filters themselves is shown
(from Fred Harris presentation)
The derivation of the band-edge filters for the FLL block in GNU Radio has based on an RRC match filter response. In the flowgraph below, the FLL block can perform course carrier recovery for all NRZ formats (L/M/S). However, the block fails to do so when the data is encoded in Manchester/Biphase line code. The Manchester line coder/decoder have been verified to be functionally okay so encoding issues can be ruled out. To add more confusion, the constellation diagram is quite okay even though no output frames are received. Could the nature of the Manchester pulse (shown below) be the reason for this problem? Should I re-derive the band-edge filter responses based on the Manchester pulse[Though to my understanding the match filter response is still the same regardless of the line code being used]?

Regards,

M.


texts: consecutive filtered blocks with filtfilt

texts:
I'm doing basic signal processing in Python. I'm filtering consecutive blocks of a signal with scipy's filtfilt (to get zero phase). I would like to concatenate several blocks after filtering, however, each filtered block doesn't seem to be properly continuous with the other filtered blocks. Is there some way to work this out?

I'd like to note that when I'm doing the same with scipy's lfilter with "Initial Conditions"(zi) parameter being passed to it at each iteration the filtered blocks seems to be continuous.

Any help will be highly appreciated.

Mike


texts:
LoRa integrate data into symbols by selecting the starting frequency of an up chirp. The resulting up chirp traverses across BW and arrives at the same place it started. At the demodulator, this chirp is multiplied by a down chirp followed by the fft to extract the data bits.
I tried to do this myself on paper but I can't help but find two different fft bins. I understand if both -ve & +ve bins are same, the demodulation is successful however, they don't aren't. I illustrate this at my best below.
As per DSP principles, the multiplication of received and locally generated base band chirps is the addition of the two instantaneous frequencies of the chirps in red and purple. Then, I have drawn a line which represents their summation in yellow. As we can see, the multiplication produces two products (in yellow), one positive and one negative frequency component that are not equal.

But according to literature , multiplication with a down chirp should produce only a single frequency bin, meaning, the positive and negative frequency products resulting from the multiplication must be of same magnitude. But, I see two different products here. How is that possible?
Would be so much helpful to me if someone can point out where I am wrong.


texts: Using ROC to find stability of system in specific example

texts:
I've started learning about finding the ROC from the transfer function, but I'm confused about an example.

$$H(Z) = \frac{2Z + 1}{Z^2 + Z - \frac{5}{16}}$$

I understand the poles lie at $z = \frac{-5}{4}$and $z = \frac{1}{4}$, and I assumed the system was non-causal because even though it has more poles than zeros, the unit circle does not lie to the exterior of the outermost pole(if I'm right), but because the ROC includes the unit circle, the system was stable. However, the solution to the problem states the system is both unstable and non-casual so I was hoping someone could explain this to me? I've only seen examples giving the ROC when both poles have been of the form $(z-a)$ however here one is $(z + b)$, so I thought the ROC was $\frac{5}{4} > |Z| > \frac{1}{4}$? I imagine this is where I have went wrong.

Thanks in advance.


texts: How to decide filter order in Linear Prediction Coefficients ( LPC ) while calculating formant frequency features?

texts:
I am new in signal processing and trying to calculate formant frequency features for different .wav files.

For calcuating formant frequency, I need three parameters values :

Linear Prediction Coefficients ( LPC )
root
angle

I am trying to calculate Linear Prediction Coefficients ( LPC ) using librosa.core.lpc in python. it takes two parameters:

librosa.core.lpc(y, order)

I have Y but I don't know how to calculate order, I have many .wav files and I have to set order to extract features from all the files. How to determine the order for all wav files to calculate LPC?

Next two things root and angle can be calculated easily like this :

 rts = numpy.roots(A)
 rts = [r for r in rts if numpy.imag(r) >= 0]

 angz = numpy.arctan2(numpy.imag(rts), numpy.real(rts))

 # Get frequencies.
 Fs = spf.getframerate()
 frqs = sorted(angz * (Fs / (2 * math.pi)))

But I am already stuck in the first step.

Thank you in advance!


texts:
Assuming I have samples of a colored noise.
I know the empirical Variance of those samples.
I also know the bandwidth of the filter which colorized the noise (Only the bandwidth, not the actual filter).

Now I want to estimate the variance of the noise if I filter those samples by a (Different) given filter.

Is there a way to estimate or bound from above and below the variance of the noise?

Thank You.


texts: Shannon-Nyquist theorem reconstruct 1Hz sine wave from 2 samples

texts:
Lets say I want to set the minimal sampling rate to reconstruct a 1Hz sine wave, according to the Nyquist-Shannon theorem that states that the maximum recoverable frequency is Fs/2 i.e. we must sample the signal 2 times the maximum frequency.

It seems obvious that the limit is a sampling frequency of 2Hz. Lets even say I interpret it as just more than 2Fs, lets say 3Fs. This would imply that 3 samples is enough.

So, in Matlab i generate:

f=1;%my 1Hz freq
Fs=2*f+1;%=3 =>a bit more than the Nyquist freq
t = linspace(0,1,Fs);
y=sin(2*pi*f*t);
subplot(1,2,1);
plot(t,y);
title('sin(2*pi*1*t)');
subplot(1,2,2);
plot(abs(fft(y)));
title('abs(fft(y))');

The ''sine'' does not even go back up to 0, not even mentioning that even if it did it would be a saw-tooth wave rather than a sine but i guess that is not the problem.

What am I missing, why do I need at least 4 samples rather than 2?

I think this is important to understand the theorem 'in practice'.

Although even with 4 samples the spike in the FT is wrong, it is at 3 rather than 2 (2 is the 1Hz since the first would be the DC freq)

In fact I need 5 samples to get finally this spike at 2 in the FFT! Why ?


texts: Calculating the extent of spectral leakage

texts:
As a follow up to this question -

Is there a way that I can calculate the amount of spectral leakage in a given bin, i.e., given an input amplitude, what percentage of the amplitude is present in the adjacent bins?

Or given an input frequency and the frequency range of each bin, can I calculate the amplitude present in each bin?


texts: Digital Anti Aliasing Filter For Waveform Band Limiting

texts:
I am working on some waveform generation software and I am trying to figure out what the best type of filter would be to use a an Anti Aliasing filter for my algorithm.
I am generating my waves in a "Raw" mathematical way, meaning that i am creating a ramp for a saw wave and my square wave consists of pure 1's and 0's. The issue is that this technique is inherently not band limited. I would like to implement low pass filter into the algorithm in order to filter of the harmonic being generated that are higher than the nyquist limit.

What would be the ideal filter design for this? What kind of roll of would i need?

I am doing this in c++ but i don't neccesarily need an answer in c++. just an equation or some resources would be good as well.

Thanks

edit:
I am not changing the sampling rate at any point in the chain. But i am still experiencing some aliasing somehow.

here is my code:

double _DSP::Saw_Wave(double* _frequency,_DSP::Saw_Data* _data){
    double _val = _data->_phasor;
    _data->_phasor += 2.0*(1.0/(_DSP::Sample_Rate::Samples_Per_Second()/ *_frequency));
    if (_data->_phasor>_data->_tolerance) {
        _data->_phasor-=2.0;
    }
    return _val ;
}

typedef struct Saw_Data{
        double _phasor = 0.0;
        double _tolerance = 1.0;
    }Saw_Data;

where _DSP::SampleRate::Samples_Per_Second() is the sampling rate.

the code generates the proper wave and frequency but has unwanted low frequency information when played at higher frequencies.


texts: What are the practical constraints on designing Sensing matrix in compressed Sensing?

texts:
In a typical compressed sensing scenario,  $y=Ax$, where $x$ is a sparse signal and $A$ is the sensing matrix.
To reconstruct the sparse signal $x$ from $y$,  $A$ should posses the Restricted Isometry Property (RIP). Therefore, if a matrix does not possess RIP the reconstruction of $x$ will not be possible. How can we check a matrix for RIP ?
Secondly what are the practical constraints on designing the sensing matrix $A$ ? Like should the matrix consists of real values or complex? What distribution to use for generating the matrix $A$ ? etc,.


texts: 16x interpolation using polyphase FIR - sum of the coefficients

texts:
I have an issue with polyphase implementation during interpolation. Lets assume I have 256 taps long FIR lowpass sinc filter. In order to interpolate by a factor of 2 I do use two FIR filters (polyphase filter) out of the main one since I don't need to calculate zero stuffed samples. The same goes with interpolation by a factor of 4 in which case four sub-filters (polyphase filters) are used.

However, there is something which bothers me a lot and I don't know what to think about it - the sum of the coefficients within sub-filters should always equal to 1 (if the gain is 1). That is the case with polyphase filters for interpolation factors of 2 and 4. However, if I set interpolation to a factor of 8 or greater sum of coefficients within each polyphase does no longer create a correct sum (which should be close to 1).


texts: Implementing Gaussian random variable by using a uniform random variable

texts:
I'm trying to write a C++ function that will return Gaussian random values, given their means and variances.

There is a library function rand(), which returns random numbers between 0 and RAND_MAX. RAND_MAX does not have a fixed value, but it is guarantied that it will be at least $2^{15}-1$. Its PDF is uniform.

I'm using Central Limit Theorem to transform this rand()into a Gaussian variable. What I'm exactly doing is a to call rand() for a user specified times, then add up their return values, then shift its mean to the user specified mean.

In the plotting above, I called my Gaussian random generator for $10^7$ times, and plotted frequencies of its return values. As you see, its variance is huge, since it is created by sum of a lot of other random values.

It successfully returns a Gaussian variable with a Gaussian PDF and with the specified mean value. However, the problem is its variance. I'm stuck at this point, because I don't know how to change its variance to the user specified value.

This is my code (incomplete for now; the parameter "Variance" is ignored):

template <class T>
T Random::GetGaussian(T Mean /*= 0*/, T Variance /*= 1*/)
{
    T MeanOfSum = NUM_GAUSSIAN_SUMS / static_cast<T>(2);
    T Rand = 0;
    for (uint64_t i=0; i<NUM_GAUSSIAN_SUMS; i++)
    {
        Rand += static_cast<T>(rand()) / RAND_MAX;
    }
    return Rand - (MeanOfSum - Mean);
}

Assume that NUM_GAUSSIAN_SUMS is 100, and RAND_MAX is 32767.

I want to change variance of the random variable according to the parameter of the function. My question is, how can I change variance of this random variable? How can I do it?


texts: Computation of the Inverse DCT (IDCT) Using DCT Or IFFT

texts: Is there a way to compute the inverse discrete cosine transform (type-2) by leveraging either a DCT, FFT, or IFFT algorithm? I have seen ways to compute the DCT using FFTs, and I've seen ways to compute IFFT using FFT. I can't quite find a simple example with description.


texts: How to generalise the Fourier transform?

texts:
The Fourier transform takes a signal and splits it into a series of sine and cosine waves.

I am told that it's supposed to be possible to split a signal into some other set of functions. My question is: How do you do this?

I'm presuming that the set of functions you use would have to have certain properties for this to work. (E.g., you have to have "enough" different functions to capture all the information of the original signal.) How do you figure out whether your set of functions is suitable? And then how do you do the actual splitting?


texts: Transfer function estimation from logarithmic sweep frequency response data

texts:
I have complex frequency response data (of an analog system) in the range of 100 Hz to 100 GHz, and it is sampled in frequency with logarithmic spacing. I would like to be able to turn this into a filter in MATLAB such that I can multiply it with the fft of time-domain signal.

I'm not seeing a good match when using invfreqs.  What is the best way to do this?


texts: Estimating the Impulse Response of the Room Using Sweep Signal Microphone Recorded Signal (Input & Output of a Convolution)

texts:
I played this signal A (a 20Hz to 20000Hz sinusoidal sweep in 10 seconds) with a studio monitor speaker in a big church, and I recorded the result B with good microphones.

The result is very reverb-ish, that's exactly what I wanted to catch.

Now a software (such as Deconvolver but non open-source) can build an impulse response from A + B, that can be later used in a convolution reverb.

It works well. But I would like to learn how to do this myself via DSP / programming.

How can I use a sweep (signal A) + recorded output (signal B) to get an impulse reponse?

Edit: In other words, if a is the original sweep, b the recorded output, and h the impulse response, how to get h from

$$a * h = b$$

Is this formulation correct? is the solution $h=a^{-1} * b$, where $a^{-1}$ is the inverse of $a$ for the convolution? How to compute a convolution-inverse of a discrete signal?


texts: Importance of phase of time-series for PSD

texts:
I have $m$ complex-valued time-series $z_m(t_n)$, which I am currently analyzing.

I started by calculating their PSDs and CSDs. These exhibit an unexpected peak at a given frequency, which is present for all PSDs and CSDs, signaling some kind of coherence between the time-series.

Interestingly this peak disappears, when I calculate the PSDs or CSDs of only the modulus $|z_m(t_n)|$ of the time-series, suggesting that the phase $\phi[z_m(t_n)]$ of the time-series is somehow related to the time-domain and thus is important for determining its the power.

My question is acutally two-fold:

1) What does the relevance of the phase of $z_m(t_n)$ vs. $|z_m(t_n)|$ tell me about $z_m(t_n)$? I.e. is there some mathematical relevance/relation that could explain this?

2) Could this be due to some kind of "principal phase" of the $z_m(t_n)$ that, by throwing it out when calculating the modulus, "decorrelates" the spectral information of my time-series in the frequency-domain?

I am a little puzzled by the behavior and hope somebody could suggest, where to look or if there is some way to unravel the relevance of the phase in this.


texts:
For linear equalizers, the weights of the equalizer are set to achieve minimum error between the equalized signal and the original input signal. First, a training sequence is sent as input, and the weights that achieve the minimum error between the equalized and the original training sequence are obtained. My question is:

In case of non-adaptive equalizers, are the weights obtained for the training sequence used for any other input signals?
In case of adaptive equalizers, how the weights are updated according to the new input signals? The error should be minimized between the equalized signal and what? There is no known sequence in this case.


texts: Energy normalization across wavelet subbands

texts:
According to the following reference: A Really Friendly Guide to Wavelets, © C. Valens, 1999
Equation 3, the wavelets are generated from the mother wavelet by scaling and translation. S is the scaling factor. The equation is divided by the square root of S for energy normalization across the different bands. I can not understand why this division achieves energy normalization.

According to the same reference, equation 14,  From Fourier theory we know that compression in time is equivalent to stretching the spectrum and shifting it upwards. I can understand spectrum stretching in this equation, but where is the spectrum shifting in the equation?1


texts: During MFCC process FFT is done to get the power spectrum but last step is DCT in order to get back to time domain. Cannot DCT be used in the first case to convert the signal to frequency domain?


texts: Masking Threshold in psychoacoustics model

texts:
I am trying to calculate masking threshold using power spectral density but is is given in the attached document that global masking threshold is the sum of individual masking threshold and threshold in quiet.Why we are adding threshold in quiet??

http://103.71.168.43:81/mis/files/coursefiles/AY2018%20-%202019/ECE/AP5074/cf~Masking~69~AP50741542606821-1.pdf


texts:
I want to convert a .wav file with a sampling frequency of 44100 Hz to a 16 depth .pcm.

I don't why I'm getting those peaks at the beginning of  .pcm plot (second figure below 2). If you could explain and tell me how I can correct it, I would very much appreciate it.

This is my attempt:

% determine the least common multiple (lcm) of fsin and fsout
fsin  = fs;
fsout = 22050;
m     = lcm(fsin, fsout);
% determine the up and down sampling rates
up   = m/fsin;
down = m/fsout;
% resample the input using the computed up/down rates
x_22 = resample(x, up, down);
audiowrite([a_filename(1:11),'_22050','.wav'], x_22, fsout);
precision = 'int16';
fidr      = fopen([a_filename(1:11), '_22050','.wav'], 'r');  % open .wav file to read
fidw      = fopen([a_filename(1:11), '_22050','.pcm'], 'wb'); % open .pcm file to write
w         = int16(fread(fidr, inf, precision));% read wav file
fwrite(fidw, w, precision);
fclose(fidr);
fclose(fidw);

fidr2     = fopen([a_filename(1:11),'_22050','.pcm'], 'r');
[data, ~] = fread(fidr2, 'short');
fclose(fidr2);

% plot
figure(1)
set(1, 'color', 'w')
subplot(311),plot(x)
grid on, box on, axis tight, title('.wav (f_s = 44100 Hz)')
subplot(312),plot(x_22)
grid on, box on, axis tight, title('.wav (f_s = 22050 Hz)')
subplot(313), plot(data)
grid on, box on,  axis tight, title('.pcm (16 bit depth)')

  These are the resuts:
  Here are all the steps of the conversion:

This is the zoom in on the thrid subplot to see the peaks:


texts: Is it numerically more stable to implement filtering as multiplication or convolution?

texts:
I'm writing a program to filter a 20,000-sample signal with a fifth-order Butterworth filter offline. I've got access to an FFT implementation. There seems to be two alternatives for implementing the filtering:

convolving the signal with the impulse response in the time domain, or
multiplying the signal with the impulse response in the frequency domain and inverse-transforming the result

These methods would be identical in the theoretical FT case. Doing it in real life with the DFT, though, I suppose things are different. Is one of the methods numerically more stable? Are there any other problems that I should be aware of? The number of calculations is not important.


texts: How to conserve symmetry of discrete spectrum after applying fractional delay in the frequency-domain?

texts:
Imagine a real-valued, discrete-time signal $h[n]$ with $N$ samples and a sampling period of $T_s = 1/f_s$. We know that the DFT of the signal $H[k]$ is conjugate symmetric because the signal is real-valued:

$H[N-k] = H^*[k], k=1,..,N-1$ and $H[0] \in R$

Now you want to apply a fractional delay $\tau [s]$ (i.e., any real-valued delay, not limited to an integer times the sampling period) in the frequency-domain:

$H_{\tau}[k] = e^{-\frac{j2 \pi k}{N}f_s\tau}H[k]$

Now verify symmetry again:

$H_{\tau}[N-k] = e^{-\frac{j2 \pi (N-k)}{N}f_s \tau}H[k] = e^{-j2 \pi f_s \tau}e^{+\frac{j2 \pi k}{N}f_s\tau}H[k]$

which means that the symmetry will only hold for $f_s \tau \in Z$ (i.e., the delay should be an integer times the sampling period). This also implies that the inverse DFT of my delayed spectrum will not be real-valued anymore, which actually makes no sense.

Therefore, my question is whether there's a straightforward way to apply a fractional delay in the frequency domain.

One way I can come up with is to apply the delay on one side of the spectrum, and then impose symmetry; however, this does not seem to yield the desired signal.

Another way is to multiply with the DFT of a fractionally delayed sinc-function, but then I am actually using a time-domain technique. Because the DFT of a sinc-function is actually a rectangular window, and delaying this sinc-function results in a phase shift in the frequency domain as described above, I understand that this way is some small correction to the first-described method which resulted in a complex-valued signal. However, I do not see what this small correction actually is in the frequency-domain, such that this technique does end up with a real-valued signal.

A third way I came up with is to do a DTFT instead of a DFT, apply the delay, and then sample this DTFT to end up with a DFT. However, I do not see how this works numerically.


texts: SNR contains $\log N$ term where $N$ stands for number of samples

texts:
Sampling the signal$ $N times increases the signal energy by a factor of $N^2$ and the noise energy by a factor of $N$. Why? This explanation is written for
              SNR(dB) = signal peak(dB) – noise floor(dB)- $10\log N$


texts: How to do smoothing without affecting phase

texts:
I have signals captured from two channels of a measurement system. There is phase and magnitude differences between these two signals.

I want to apply nonlinear smoothing to the magnitude, but the phases should remain the same after the process.

This is what I have done both channels:

Perform FFT
calculate magnitude
calculate phase
apply nonlinear smoothing to the magnitude
get real and imaginary part from smoothed magnitude and phase
perform inverse FFT.

However the resultant signal is not as desired. It appears as a two sided time signal or has some additional signals in the end.

Is there a solution to this problem?


texts: Computation of parameter filter to match a given frequency response

texts:
I'm looking for the practical way to compute the parameters ($a_1$, $a_2$) of a digital filter to match a certain frequency response. I'm studying the 12-poles filter of the vintage component SP0256 which is used in speech synthesis, using bandpass filters to create frequency response of allophones.
The 12-poles filter is composed of six 2-poles filters :

$$H(z)=\frac{1}{1-2 a_1 z^{-1}-a_2 z^{-2}}$$

I came accros a github of the FPGA implementation of the SP06256 which provide a matlab example of this filter for the 'EH' allophone. I ploted the 6 frequency responses in blue and the final one in red with matlab :

I could find the formant frequency for the other allophones, but I can't figure out how to compute the filter parameters to match the desired frequency response.


texts:
I have a design in which a signal from ADC is filtered and then played by DAC. In FPGA the real signal is converted to quadrature using DDS with 70 MHz frequency. After filtering, the baseband signal is sent to DAC. DDS of DAC is also set to 70 MHz.

The problem is that there is always a spur in output of DAC at 70 MHz. Noise level is almost -93 dBm and this spur is -82 dBm while the main signal level is -15 dBm.

Level of spur does not depend on the signal level, as long as signal passes through filter. My first guess was that this is because of signal truncation. There are two steps of truncation in design. I have used convergent rounding instead of truncation but the spur is still there; decreased by almost 6 dB but still there. then I wondered that this might come from ADC, but changing the frequency of input signal does not change the spur; It's always at 70 MHz. Then I taught that it might be because of DAC LO but when the signal is in rejection band of filter, the Spur disappears.

Do you have any idea what makes this spur and how can I remove it?

Any help would be appreciated.


texts: Why do I have to divide by $\pi$ in order to scale amplitude of sine frequency spectrum?

texts:
Recently I found my old code to transform sine wave from time domain to frequency domain. Below I post simplified version of it.

from math import sin, pi
from matplotlib.pyplot import *
from pylab import zeros, fft

fs = 500
duration = 2 * pi
amplitude = 2
frequency = 25
values_in_time_domain = []

for t in range(int(duration * fs)):
    values_in_time_domain.append(amplitude * sin(frequency * t / fs))

subplot(2, 1, 1)
title('Sinus from 0 to 2 pi in time domain; A = 2, f = 25, fs = 500')
plot(values_in_time_domain)

transform = fft(values_in_time_domain)
transform_abs = [abs(x) for x in transform]
values_in_freq_domain = [transform_abs[f] / fs / pi for f in range(int(fs / 2)+1)]

subplot(2, 1, 2)
title('Above sinus in frequency domain')
plot(values_in_freq_domain)

show()

Why do I have to divide every element from transform_abs by pi in the following: transform_abs[f] / fs / pi?
I noticed I need to do it in order to have properly scaled $O_y$ axis.


texts: Solve equations with conjugates multiplications

texts:
If I have a two variables $x_1$, $x_2$ and two equations. In the first equation the first variable $x_1$ is multiplied with the conjugate of same number which is multiplied with $x_2$, and in the second equation, the exact inverse happens.  For example, we have a number $1-2j$ so:
$(1-2j)x_1 + (1+2j)x_2 = y_1$  ... (1)
$(1+2j)x_1 + (1-2j)x_2 = y_2$  ... (2)
I'm confused, what's the proprieties of two above equations, I mean can we solve $x_2$ and $x_1$ if  $y_2$ and $y_2$ are known? .. If not, can we get the value of $y_2$ based on the equation (1), it means we still have two variables but with one equation, we need to get the second equation which is very similar to equation (1)


texts: Opencv calcHist and calcBackProject in plain english

texts:
i've started experimenting with c++ and opencv because of i want to learn image processing.

Now, my first exercise is to create a skin detector with calcHist and calcBackProject.

But i don't understand few things:

the statistic interpretation of back project, and why it is named "back project"
i have quite good comprehension of what ranges parameter is in calcBackProject. But i'm really stucked with ranges parameter in calcHist function.
For more detection precision i think could be a good thing use backgrojection in more levels: in each channel of r-g-b ans in each channel of h-s-v. But i don't know how i can combine the different results of calcBackProject of separate channels of rgb and hsv.

And i think that my not-well comprehension is caused by lack in theory of what i'm doing with those 2 methods (see the first point). So please explain me in plain english.


texts: Confusion on using FFTW 3.3.4 (C version) for real and even data sets

texts:
I am using FFTW3 to transform the autocorrelation function of a particular function. I expect the transform to be real since the autocorrelation function is symmetric with respect to negative and positive lags and its Fourier transform give the power spectrum. However, I get nonzero complex part and negative real parts after the transformation. What I may be doing wrong?

I am providing the code below:

/* Compilation and execution command:
 * gcc q1part2.c nrutil.c qtrap.c trapzd.c -lfftw3 -lm -o q1part2 && ./q1part2
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <fftw3.h>
#include "nr.h"
#include "nrutil.h"
double h = 0.01;
double a = 0.0, b = 10.0, p;

double func(double t)
{
    return  cos(10 * t) + sin( 20 * t);
}
double intgrd(double t)
{
    return func(t) * func(t + p);
}

int main()
{
    int  N = round(((b - a) / h) + 1);
    double y[2*N - 1], t[2*N - 1]; /* t array hold time delays */
    int i;
    fftw_complex *in, *out;
    fftw_plan MYPLAN;
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
    MYPLAN = fftw_plan_dft_1d(2*N - 1, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    for(i = 0; i < 2*N - 1; i++)
    {
        /* The definition of autocorrelation from the book "A First Course on
        * Computational Physics" by Paul DeVries was used.
        */
        p = i * h - 10.0;
        t[i] = p;
        if(p >= 0)
        {
            y[i] = sqrt(1 / (2 * M_PI)) * qtrap(intgrd, a, b - p);
            in[i] = y[i];
        }
        else
        {
            y[i] = sqrt(1 / (2 * M_PI)) * qtrap(intgrd, -p, b);
            in[i] = y[i];
        }
        //printf("%f \t %f \n", t[i], y[i]);
    }
    fftw_execute(MYPLAN);
    for(i = 0; i < N; i++)
    {
        printf("%f \t %f  %f \n", i * 2*  M_PI / (2 * (b - a)), creal(out[i]),cimag(out[i]));
    }
    fftw_destroy_plan(MYPLAN);
    fftw_free(in);
    fftw_free(out);
    return 0;
}

The function qtrap is a numerical integration algorithm provided in the code from the second edition (plain C) of the book Numerical Recipes. I am really confused about the usage of the library and would like to get pointed in the right direction.

Thanks in advance


texts: Highpass filter design that is better than the Moving Difference filter

texts:
The question is actually from the problem set 1.3 of the "Wavelets and Filter Banks" book (authors Nguyen/Strang).
Problem 4: Invent a Highpass filter K with 3 or 4 taps (coefficients) that is better than the moving difference H1. The goal is:
$$
 \left| K(w) \right| < \left| H_1(w) \right| \;\;\; for \;\;\; 0 < \left| w \right|  < \pi/2
$$
And,
$$
 \left| H_1(w) \right| < \left| K(w) \right| < 1 \;\;\; for \;\;\; \pi/2 < \left| w \right|  < \pi
$$
We know that
$$
H_1(ω) = \frac{1}{2} (1-e^{(-jω)} ) \\
\left| H_1(ω) \right| = \left|sin(ω/2)\right|
$$
I'm not sure if the best way to tackle this is by going for an anti-symmetrical filter (4 taps, k1 = -k3, and k2 = -k4). Should I just assume coefficients until I actually meet the condition 1 stated? Or is there another approach? Perhaps setting an equation system to find out about h1 and h2 for K(w)?
Also, does condition 2 mean that for this interval of |W|, the filter needs to be invertible too?
Any help is dearly appreciated.


texts: RMS calculation in frequency domain after windowing

texts:
I can calculate RMS in frequency domain as derived from Parseval's Theorem.

But what if I have applied a windowing function before doing the FFT (in my case a Hann window)?. Now the RMS values are not correct any more. How do I calculate the correct RMS values?


texts: What is the real meaning of "frequency" for a non periodic signal?

texts: For a non periodic signal, in the first one second  the frequency of a particular amplitude (value) may be of 4 Hz, and then in the next one second  it may or may not be 4 Hz. Then which value of frequency should you consider?


texts: Comparing spectra from hyperspectral images

texts:
I extracted the mean spectra of regions of interests of different tissues from  hyperspectral images. How to compare those spectra by the difference of their reflectance and their spectral behavior  (I mean the spectra look similar but are shifted to the left or right).

Example of the mean spectra of two regions:

I would like to know the similarity between the two spectra. I tried cross correlation and Root mean square error but the results were not satisfying.

The goal is to predict a genetic distance of species from the spectra. For instance if the genetic distance of species 1 and species 2 is very height I would expect a different spectra. If the distance is 0 or very low,m the spectra are the same /very similar.


texts:
I am computing the link budget for a communication-radar (CoRadar) system. The formula for maximum distance, in communication, is:

$d_{max} = \sqrt{\frac{P_t}{\tau}\frac{G_t G_{r, comm} \lambda^2}{(4\pi^2)FkT_0 B_w SNR_{comm}}}$

where,

$P_t$ is the transmitted power

$\tau$ is the duty cycle

$G_t$ is the transmission antenna gain

$G_{r, comm}$ is the communication receiver antenna gain

$\lambda$ is the wavelength

$F$ is the noise factor

$k$ is the Boltzmann's constant

$T_0$ is the noise reference temperature

$B_w$ is the bandwidth

$SNR_{comm}$ is the $SNR$ of the communication system

I am using a LFM chirp signal that sweeps from $3.5\times10^6$ to $8.5\times10^6$. Which would be the frequency, which would be associated to a wavalength $\lambda$, that should I use for the formula? $8.5\times10^6$? Or, $\frac{8.5\times10^6 + 3.5\times10^6}{2} = 6\times10^6$ (i.e. a middle point)?

And which would be the bandwidth? $8.5\times10^6 - 3.5\times10^6 = 5\times10^6$?

Many thanks.


texts: FFT Matlab - Meaning of Frequency Vector

texts:
I'm following a tutorial about the FFT. It's well explained but I don't understand the meaning of the frequency vector:

Fs = 150; % Sampling frequency
t = 0:1/Fs:1; % Time vector of 1 second
f = 5; % Create a sine wave of f Hz.
x = sin(2*pi*t*f);
nfft = 1024; % Length of FFT

% Take fft, padding with zeros so that length(X)
is equal to nfft
X = fft(x,nfft);
% FFT is symmetric, throw away second half
X = X(1:nfft/2);
% Take the magnitude of fft of x
mx = abs(X);
% Frequency vector
f = (0:nfft/2-1)*Fs/nfft;

% Generate the plot, title and labels.
figure(1);
plot(t,x);
title('Sine Wave Signal');
xlabel('Time (s)');
ylabel('Amplitude');
figure(2);
plot(f,mx);
title('Power Spectrum of a Sine Wave');
xlabel('Frequency (Hz)');
ylabel('Power');

Can someone explain me what is this for:

f= (0:nfft/2-1)*Fs/nfft;

Here is the link


texts: Recording and playback of sound, frequencies question

texts:
If I record a sound at a rate of 8000 samples per second I know that the highest frequency that can be reproduced from the sound is 4000 Hz. What I want to know is;

  If the sound contains frequencies above 4000 Hz when I play the sound
  back what would I hear and what frequencies would be produced at the
  points where the frequency of the original sound is greater than
  4000Hz.


texts: Do FIR filter have any application in control theory?

texts:
I know about IIR filter as discrete pendent to transfer functions in the Laplace domain. So it is actually quite simple to convert the function of the control system and finally receive the discrete output vector. (Bilinear Transform...)

But, how about FIR filter? Can we use them  to approximate functions of control systems.


texts: Finding similar fragment in recorded speech

texts:
I use MFCC algorithm for calculation of Mel-frequency coefficient for 2 files and DTW algorithm with Euclidean distance for finding similarity. I have results values around 400-600 on the files with the same word (cut from recorded speech of the same person). Is it right? Or I have a mistake in my implementation? I didn't find an implementation of MFCC in C#, so I edited some files from comirva (java), can it be the problem?
I have also tried Chi-Square distance, but it gives a negative results.


texts:
Let's suppose i want to implement a Buttersworth filter of order 12 (N=12).
First of all ,I can decompose the Buttersworth polynomial in elemental polynomials:

Now, if i know how to physically implement the elemental polynomials, how can i construct the original Buttersworth filter of order 12?


texts:
Suppose I have an IIR filter in the $z$-domain in the following form:
$$
H\left(z\right)=\frac{1}{1-0.2z^{-1}-0.1z^{-2}}
$$
How do I represent this in MATLAB?
I am pretty sure if I just listed the coefficients as H = [1 -0.2 -0.1];
this would be considered wrong. Can somebody help me with this please?


texts: why bandpass filtering (butterworth) centers the signal

texts: I run bandpass filtering via the butter() in matlab to apply butterworth bandpass filtering to find out the signals are centered at zero after the filter. Would this always happen for any kind of low/high-pass filtering and why?


texts:
I am new to signal processing and I want to calculate the 3 dB - cut off frequency of this FIR filter, which has the following transfer function:

For the evaluation i need to replace the z with  and set the squared magnitude equal to 1/2.

But the problem is when I want to solve this for omega zero is the result. My question is what am I doing wrong here, or is something missing ?


texts: Relating bandwidth and wavelength of a carrier

texts:
This came from a homework question, which I've solved. Just need some clarification.

I am asked to determine the bandwidth in a given amount of spectrum at a given wavelength. To find this, I know that;

$$\delta_f = \frac{c\cdot \delta_\lambda}{\lambda^2}$$

where $\delta_f$ is the bandwidth, $\delta_\lambda$ is the amount of spectrum, $\lambda$ is the given carrier wavelength, and $c$ the speed of light.

That's all well and good, but I can not find any reason for why this formula works, or where it comes from. Where does this formula come from, or what is it called, or how is it derived?

The question reads as follows;

  How much bandwidth is there in $0.1 \textrm{ micron}$ of spectrum at a wavelength of $1 \textrm{ micron}$?


texts:
I'm applying a radial low-pass filter on a $H{\times}W$ pixel 2D image $I$ where $H \neq W$. To implement this filter, I apply 2D DFT to $I$ to produce $I'$, the $H{\times}W$ frequency domain representation of $I$. I apply an element-wise multiplication $I' \odot M$ with a filter matrix $M \in \{0, 1\}^{H{\times}W}$. Let $u \in \{0, H-1\}, v \in \{0, W-1\}$ be indices in $M$. $u$ represents cycles/entire height of image, while $v$ represents cycles/entire width of image. An ideal low-pass filter with radius $R_{img}$ is given by
$$\sqrt{u^2 + v^2} < R_{img},$$
whose units are in cycles/image. I want to change the units of the filter parameter $R_{img}$ from cycles/image to cycles/mm, given that $\gamma$ = mm/pixel.

Here is what I have so far, please let me know if this is correct. Since the filter radius can't exceed $\min(H, W)$, let $f_s = \min(H, W)$ denote the sampling frequency of the image across both dimensions in terms of cycles/image, and $u, v$ denote row and column indices of the $H{\times}W$ representation of the image in frequency domain after applying shifted 2D DFT. We then have

$$\sqrt{u^2 + v^2} < R_{img}\\
\frac{1}{f_s \gamma} \sqrt{u^2 + v^2} < \frac{R_{img}}{f_s \gamma}\\
\sqrt{\left(\frac{u}{f_s \gamma}\right)^2 + \left(\frac{v}{f_s \gamma}\right)^2} < \frac{R_{img}}{f_s \gamma},$$

where the units have changed from cycles/image to cycles/mm. In other words,
$$R_{mm} = \frac{R_{img}}{f_s \gamma}.$$


texts: Estimate Parameters of Linear Combination of Harmonic Signals with Partial Known Frequencies

texts:
Consider real signal of the form:
$$
\begin{align*} x(n) = &\sum_{m=1}^{M} \left( c_m \cos(2\pi f_m n) + s_m \sin(2\pi f_m n) \right) =\\
&\sum_{m=1}^{M} a_m \cos(2\pi f_m n + \phi_m)
\end{align*}
$$
where signal values of $x(n)$ are known for $n = 1,\dots, N >> M$ and some frequencies $f_m$ are known real numbers in $(0.0,0.5)$.
What is the best way to estimate $c_m$ and $s_m$ real parameters?
Estimation of the form:
$$c_m + j s_m = 2 \frac{\displaystyle\sum_{n=1}^{N} x(n) \exp(2 \pi j f_m n)}{N}$$
is not very accurate ($j^2=-1$), but windowing (e.g. cosine window) improves estimation accuracy significantly:
$$c_m + j s_m = 2 \frac{\displaystyle\sum_{n=1}^{N} w(n) x(n) \exp(2 \pi j f_m n)}{\displaystyle\sum_{n=1}^{N} w(n)}$$
Why window is so effective here? Is it some kind of well-known formula?
The least-squares method also gives good results if frequencies are known for all terms.
Are there any other good methods for the case when not all frequencies are known? Perhaps based on SVD or some kind of orthogonalization.
Test code in Mathematica:
F1 = 0.123000 ;
F2 = 0.312874 ;

LENGTH = 1024 ;
RANGE = Range[LENGTH] ;

SIGNAL = 0.0 ;
SIGNAL = SIGNAL + 1.0*Cos[2*Pi*F1*RANGE] + 0.5*Sin[2*Pi*F1*RANGE] ;
SIGNAL = SIGNAL + 0.05*Cos[2*Pi*F2*RANGE] + 0.01*Sin[2*Pi*F2*RANGE] ;

WINDOW = HannWindow[N[Rescale[Range[LENGTH],{1,LENGTH},{-1,1}]]]^2 ;

(* NO WINDOW *)
2*(Total[SIGNAL*Exp[2*Pi*I*F1*RANGE]]/LENGTH)
2*(Total[SIGNAL*Exp[2*Pi*I*F2*RANGE]]/LENGTH)
(* OUTPUT: 0.9994929536863852`\[VeryThinSpace]+0.49994639046951617` I *)
(* OUTPUT: 0.04819835064062154`\[VeryThinSpace]+0.012127390241918076` I *)

(* WITH WINDOW *)
2*(Total[SIGNAL*WINDOW*Exp[2*Pi*I*F1*RANGE]]/Total[WINDOW])
2*(Total[SIGNAL*WINDOW*Exp[2*Pi*I*F2*RANGE]]/Total[WINDOW])
(* OUTPUT: 0.9999999999775656`\[VeryThinSpace]+0.5000000000025541` I *)
(* OUTPUT: 0.05000000004559643`\[VeryThinSpace]+0.009999999958900167` I *)

(* LEAST SQUARES *)
(* ONE KNOWN *)
LeastSquares[{Cos[2*Pi*F1*RANGE],Sin[2*Pi*F1*RANGE]}//Transpose,SIGNAL]
(* ALL KNOWN *)
LeastSquares[{Cos[2*Pi*F1*RANGE],Sin[2*Pi*F1*RANGE],Cos[2*Pi*F2*RANGE],Sin[2*Pi*F2*RANGE]}//Transpose,SIGNAL]
(* OUTPUT: {0.999957445314386`,0.4999499273288035`} *)
(* OUTPUT: {1.`,0.49999999999999983`,0.05000000000000022`,0.010000000000000056`} *)


texts: What are the mathematical tools that I can employ for automatic music composition?

texts: I would like to do research on automatic music composition. At the lowest level , I would like to start from discrete frequency levels as states. Here I would like to know the mathematical tools I can use for this? I have heard of markov chains and chaos theory etc used in this purpose. But one professor recently told that graph theory will be helpful in automated music composition. Is it so? Can someone recommend me the books required for starting such a research.


texts:
In an OFDM receiver, how could we find the start of the symbol for FFT operation.

I know the CP-based method to find the start of symbol. However, what if the sampling is off by some duration. We will not get the same digital signal at the transmitter. How serious is that issue? Is it simply addressed in the channel estimation in the frequency domain?

Thanks in advance.


texts: Fixed point channel model for OFDM system

texts:
I'm testing an OFDM system implemented in fixed point. The data format is Q11. My system work's fine but I need to test it under some channel for evaluation of the design before field testing. The system will be tested on Rician fading channel but for now I'm working on AWGN testing first. I asked a similar question here:
Apply AWGN noise to QPSK-OFDM symbol

But I realized that the problem is in the fixed point - float point conversion. I had a simulation of my system using 32b and the results agree with the theoretical BERXEbN0 for QPSK systems.

I'm using the following code for apply noise:

def awgn(self,data,noise_power):
    data_dq = data*(2.**(DATA_WIDTH-1))
    Es = sum(abs(data_dq)**2.)
    Eb = Es/(self.OFDM_SIZE*2.)
    EbN0_dB = noise_power + 10*log10(self.t_symbol/(self.t_symbol - self.t_cp))
    EbN0 = (10.)**(EbN0_dB/10.)
    N0 = Eb/EbN0
    noise = sqrt(N0/2.)*(random.standard_normal((data_dq.size,)) + random.standard_normal((data_dq.size,))*1j)
    return around(data_dq + noise)/(2.**(DATA_WIDTH-1))

where DATA_WIDTH = 12.

With this implementation the system agree for some regions of EbN0 curve and goes far from the curve in the other regions. With this number of bits I'm getting $10^{-2}$ BER with an EbN0 of 10dB.

My question is: how to handle the conversion between fixed and float point in this context?

update

I had changed some parts of the above code:

def awgn(self,data,noise_power):
    data_dq = data
    Es = sum(abs(data_dq)**2.)
    Eb = Es/(self.OFDM_SIZE*2.)
    EbN0_dB = noise_power + 10*log10(self.t_symbol/(self.t_symbol - self.t_cp))
    EbN0 = (10.)**(EbN0_dB/10.)
    N0 = Eb/EbN0
    noise = sqrt(N0/2.)*(random.standard_normal((data_dq.size,)) + random.standard_normal((data_dq.size,))*1j)
    return data_dq + around(noise*(2.**(DATA_WIDTH-1)))/(2.**(DATA_WIDTH-1))

I guess it will not affect the system performance but it's here for keep updated.

Below are the curve from my last simulation:

As you can see I have a floor in my present design.

I'm asking if this floor is due to my fixed point system performance or it's an effect of the channel model.

Update[Solution]

Just fro add an update to the problem for future reference. The problem was in the parameters of the fixed-point FFT core used. The effects in the figure above are due to clipping, as stated from the helpful comments in this topic. Thank you for your support.


texts: I have 24-bit data in two's complement 32-bit integers. I need Q31 data

texts:
I have been trying to learn about Q31 format (as well as Q15), but I am having trouble understanding how it works.

So my CODEC is giving me 24-bit two's complement data. I have converted this data successfully to 32-bit two's complement data and back. I would like to so some math on this data using Q31 format. For starters, hopefully something like attenuation. I am using a PIC32MZ series processor. Here is some of the code I am working with.

#define Q31 int
#define Q63 long long

Q31 Q31multQ31 (Q31 a, Q31 b){
    Q63 y;
    y = (Q63) a * (Q63) b;
    y >> 31;
    return ((Q31) y);
}

// Super basic attenuation (hopefully)
left_output =  Q31multQ31(left_input,0b00000000001111111111111111111111);

This seems to give me awful noise. left_input and left_output are int (signed 32-bit). Is there some sort of conversion I am missing to use these as Q31?


texts: How to implement IDFT using the positive frequencies of DFT, modulation and resampling on MATLAB?

texts:
This question is a continuation of the discussion at How can a linear operator on DFT vector produce the same vector using only half of the DFT vector?, where the higher level motivation is better explained, in case there is interest.

Question:

How to take the positive frequencies from an hermitian symmetric DFT vector, take the IDFT of this (which is complex), multiply by a carrier to shift the spectrum, resample, take the real part of it, and obtain the same IDFT result as the result that would be obtained by the IDFT of the original hermitian symmetric symbol?

Details:

I would really appreciate to have my understanding and the concepts I mention checked.

Consider that the DFT vector $\mathbf{X}$ of a certain signal is hermitian symmetric and represents a bandwidth of $B$ Hz. If I consider a signal represented by the positive frequencies (from DC to Nyquist), it is the same as shifting the positive side of the spectrum, and centralizing it at frequency 0. Now, the new spectrum is not hermitian symmetric anymore, yields a complex IDFT and represents a Bandwidth of $\frac{B}{2}$.

To use a numerical example, for $N=8$, considering from DC to Nyquist, that would give 5 positive frequencies and 3 negative. Then, I would shift the 5 DFT samples (correspondent to the positive frequencies) to the center of the spectrum. In this case, the new DFT size would be odd, and the Nyquist frequency (at k= 5/2 = 2.5) wouldn't fall on a DFT sample

At this point, I have a signal bandlimited to $\frac{B}{2}$. Then, if I multiply this signal with a carrier of frequency $\frac{B}{2}$, I shift the spectrum again to the positive and negative side, such that the original spectrum is reconstructed. However, since this operation would increase the bandwidth back to $B$ Hz, the product $x[n]*c[n]$ would need to have its sampling frequency increased by a factor of 2, so that the frequency spectrum can actually represent the spectrum of the signal without aliasing.

The carrier should be an exponential:

\begin{align}
c[n] = e^{j2\pi n k_0/N}
\end{align}
where $k_0$ is the number of DFT samples to shift. In this case, $k_0 = \frac{N}{2}$, such that:
\begin{align}
c[n] = e^{j \pi n}
\end{align}

Finally, after doing resampling, I should take the real part of the signal and that should be equal to the original IDFT (the one obtained by simpling taking the IDFT of the original hermitian symmetric symbol).

However, I am confused mainly by the following issues:

From DC to Nyquist there are 5 DFT samples. That would yield a complex IFFT with size 5, Then, I would multiply point-wise by 5 samples of the carrier and upsample by 2, resulting in 10 samples. On the other hand, the original IDFT has length 8.
I'm not sure about the rationale behind taking the real part.
I'm not sure if the exponential is the correct carrier.

MATLAB implementation:

Here is the code I'm trying to make work.

N=4;
A=rand(N-1,1) + j*rand(N-1,1); % Random positive carrier
X=[0; A; 0; flipud(conj(A))]; % Hermitian symmetric DFT vector

x = ifft(X)

% Positive side of the spectrum shifted to the center
Xp = [X((N/2 + 1):N+1); X(1:N/2);];

% Complex IFFT:
xp = ifft(Xp);

% Carrier
n=(0:N);
k0 = (N+1)/2;
carrier = exp(j*2*pi*k0*n/(N+1)).';

% Modulation
x_up = xp.*carrier;

% Increase the sampling rate by a factor of two
x_re = resample(x_up,2,1);

real(x_re)

I basically want the two vectors that are printed to the workspace to be equal: x and real(x_re)


texts: What is the time-integration property in the Fourier series analysis?

texts:
In the continuous Fourier series properties for a periodic continuous-time signal, we have time-integration property:

$$
\int_{-\infty}^t x(\alpha)d\alpha \leftrightarrow \frac{a_k}{jk\omega_0}
$$

where $a_k$ is the Fourier series coefficients of the signal $x(t)$. Now I am confused. Is the above integral equivalent to $\int x(t)dt$ ?

For instance, the Fourier series coefficients of $\cos(\omega_0 t)$ are: $a_1=a_{-1}=\frac{1}{2}$, other $a_{k}=0$.

and the Fourier series coefficients of  $\sin(\omega_0 t)$ are: $b_1=\frac{1}{2j}, b_{-1}=\frac{-1}{2j}$, other $b_{k}=0$.

and we have: $\int \cos(\omega_0t)dt=\frac{1}{\omega_0}\sin(\omega_0t)$

Then:

$$
\textrm{F.S.}\left\{\int \cos(\omega_0t)dt\right\} = \frac{a_k}{jk\omega_0} = \textrm{F.S.}\left\{\frac{1}{\omega_0}\sin(\omega_0t)\right\}= \frac{1}{\omega_0}b_k
$$

Thus:

$$
b_k=\frac{a_k}{jk}
$$

as we expected.


texts: Matlab: How to implement expectation of product of hermitian matrix

texts:
I want to implement the following equation in Matlab:

$\mathbf {R}_\mathbf {z}^a(\tilde{f}) = \mathbb {E} \left[\mathbf {z}(\tilde{f}) \mathbf {z}^H(\tilde{f}+a) \right]$, $ \tilde{f} \in \left[0, f_s-a \right] $

where $ \mathbf {z}(\tilde{f}) $ is a matrix, and $ a $ takes values in the interval $ \left[0, f_s-a \right] $.
I am not sure how to implement the expectation operation, Is there some equivalent implementation in the time domain?


texts: edge detection evaluation medical test images MATLAB

texts:
I programmed an algorithm for edge detection, i want to compare it with other algorithms but with a criteria (Pratt's figure of merit, GMC1) not only a a visual comparaison.

the thing is PFOM needs an image reference (known ground truth) which i don't have, and i tried to find a GMC1 algorithm but i failed.

if any body can help me with reference images (known ground truth) (xray) or implementation would be much appreciated


texts: Why Fourier transform is not sufficient and we have to use Laplace transform?

texts:
Is there an easy way to explain the motivation behind the use of Laplace transform instead of Fourier transform?
Isn't that any periodic function can be represented by sines and cosines? - Why to introduce exponential idea?
Why not using differential equations with Fourier transform? An example would help.

*Asked the same question a while ago at math.stackexchange but no answers given.