texts: How to "scale" the FFT when using it to calculate discrete convolution?

texts:
As you probably know, the discrete convolution $ H = F \ast G $ of some $ F \left[ x \right] $ and some $ G \left[ x \right] $ can be calculated using the Fast Fourier Transform (FFT).

To do this, one would ...

Set $ M $ to the length of the sequence $ F $ and $ L $ to the length of the sequence $ G $.
Choose an $ N = 2^n $ with $ n \in \mathbb{N} $ such that $ N \ge \left(M + L \right) - 1 $.
Extend both sequences $ F $ and $ G $ to a length of $ N $ by appending zeroes (on the "right"). (This is necessary to prevent time-domain aliasing from occuring due to the fact that we perform circular convolution instead of linear convolution.)
Calculate $ \hat{F} = FFT \left( F \right), \hat{G} = FFT \left( G \right) $.
Calculate $ \hat{H} = \hat{F} \odot \hat{G} $, where $ \odot $ is the Hadamard operator.
Calculate $ H = FFT^{-1} \left( \hat{H} \right) $.

Now, the thing is that there are different variants of the DFT (and therefore FFT) and they differ from each other regarding their scaling factors.

The "standard" version is as follows.

DFT: $$ X_{k} = \sum_{n = 0}^{N - 1} x_{n} \cdot e^{\frac{-i \cdot 2 \cdot \pi \cdot k \cdot n}{N}} $$

IDFT: $$ x_{n} = \frac{1}{N} \cdot \sum_{k = 0}^{N - 1} X_{k} \cdot e^{\frac{i \cdot 2 \cdot \pi \cdot k \cdot n}{N}} $$

However, there is also an "orthonormal" version of the transforms.

DFT: $$ X_{k} = \frac{1}{\sqrt{N}} \cdot \sum_{n = 0}^{N - 1} x_{n} \cdot e^{\frac{-i \cdot 2 \cdot \pi \cdot k \cdot n}{N}} $$

IDFT: $$ x_{n} = \frac{1}{\sqrt{N}} \cdot \sum_{k = 0}^{N - 1} X_{k} \cdot e^{\frac{i \cdot 2 \cdot \pi \cdot k \cdot n}{N}} $$

And, as far as I know, in principle, one could in fact use any scaling factors to transform back and forth between time and frequency domain, as long as the product of both is $ \frac{1}{N} $.

As long as you only do one transform, then perform linear manipulations in the frequency domain, then do an inverse transform, it does not really matter, which way you scale the transforms, however, when we calculate the convolution of two vectors by taking the DFT/FFT of both, multiplying them, and then taking the inverse transform, the scaling seems to matter.

When using the "standard" transforms, both $ F $ and $ G $ get scaled by $ 1 $ (not scaled at all) on the "forward" transform, then multiplied, then get scaled by $ \frac{1}{N} $ on the inverse transform. Thus, in total, you have one scaling by multiplication with $ \frac{1}{N} $.
When using the "orthonormal" transforms, both $ F $ and $ G $ get scaled by $ \frac{1}{\sqrt{N}} $. Then the transformed functions are multiplied, so that the product of the two effectively got scaled by $ \frac{1}{N} $. Then, on the inverse transform, they get scaled by $ \frac{1}{\sqrt{N}} $ again, so that in total, you got a scaling of $ \frac{1}{N \cdot \sqrt{N}} = N^{\frac{-3}{2}} \neq \frac{1}{N} $ so that the result will be different from that obtained using the "standard" transform.

Also note that, either way, we always scale by $ N $ (or some - possibly fractional - power of $ N $), which is neither (a power of) the length of the sequence $ F $ (which is $ M $), nor (a power of) the length of the sequence $ G $ (which is $ L $) but something rather unrelated to it. In the end, we chose $ N $ so that it fits our algorithm. Most FFT algorithms require the vectors to have a size which is a power of two, so we chose one. I'm not sure if (and how) one would have to compensate for that.

Summary: How should I scale the FFTs (or apply additional scaling before/afterwards to compensate) so that $ FFT^{-1} \left( FFT \left( F \right) \odot FFT \left( G \right) \right) \stackrel{!}{=} F \ast G $?


texts: Transforming an audio signal to Mel Frequency Cepstral Coefficients is broadly used in tasks involving learning on audio. I was wondering, is this transform invertible with some good approximability properties? Are there methods to reconstruct an audio signal from its Mel Frequency Cepstral Coefficients?


texts: Multidimensionnal wavelet admissibility criteria

texts:
I've read there are two admissibility criteria for wavelets, both of which are designed to preserve total power of the signal (source: http://en.wikipedia.org/wiki/Wavelet#Mother_wavelet, as well as various scientific papers)

Condition for zero-mean: $$\int_{-\infty}^{\infty}\psi(t)dt=0 $$
Condition for square norm one: $$\int_{-\infty}^{\infty}|\psi(t)|^2dt=1,$$

where $\psi(t)$ is the wavelet kernel.
This leads to having a normalization factor of $\frac{1}{\sqrt{2\pi}}$ for the Gabor/Morlet wavelet.

Now, my question is: how do these admissibility criteria apply to 2D wavelets? (Let's say $\psi(x,y)$.)

My guess would be:

$$\int_{x=-\infty}^{\infty}\int_{y=-\infty}^{\infty}\psi(x,y)dx dy=0 $$
$$\int_{x=-\infty}^{\infty}\int_{y=-\infty}^{\infty}|\psi(x,y)|^2dxdy=1,$$

but I get inconsistent results in my computations.

Is my guess right? If not, what are they? Can you provide a reliable source?


texts: RMS value of a signal with a DC offset in MATLAB

texts:
t = 0:0.001:1;
y1 = 400 * sin(2*pi*20*t);
y2 = 200 + 400 * sin(2*pi*20*t);
plot(t,y1)
hold on
plot(t,y2,'-g')
grid on;
rms_y1 = sqrt(mean(y1.^2))
rms_y2 = sqrt(mean(y2.^2))

In the above MATLAB code there are two signals y1 and y2.

y1 is a sinusoidal signal 20 Hz with 400V amplitude.

y2 is a signal which is nothing but y1 + 200V DC offset.

So I was thinking the "rms of y2" would be equal to "rms of y1 + 200". Since rms of a constant is the constant itself.

But rms of y1 is 282.7014; and rms of y2 turns out to be 346.2948 (not 282 + 200).

What am I doing wrong here?


texts: How best to digitize a noisy phone recording for filtering?

texts: How do I digitize an analog telephone recording and filter it with a Chebyshev filter?


texts: Regularized Least Squares by Laplacian Operator - Optimal Value of the Regularization Factor (Lagrangian Multiplier)

texts:
Consider the cost function

$$f(X,\lambda) = \|AX-b\|_2^2 + \alpha \|LX\|_2^2$$

$A:$Measurement matrix($R_{m\times n}$,$m \ll n$), $b:$observation vector($R_m$), $L:$Laplacian operator($R_{n \times n}$), $X:$vector form of an image($R_n$)

(The above cost function can be read, non-mathematically as minimize

$AX=b+error$ subject to $LX=0$)

$$\frac{\partial f}{dX} = 0 \implies 2X^T(A^TA) - 2X^T(A^Tb)+\lambda(L^TL) = 0.........(1)$$

$$\frac{\partial f}{d\lambda} = 0 \implies \|LX\|_2^2=0.........(2)$$

From $(1)$ $X_s = (A^TA+\lambda L^TL)^{-1}A^Tb.........(3)$

$(3)$ in $(2) \implies \|LX_s\|^2_2 = 0..........(4)$

It is very clear from $(3),$ $\lambda$ plays a very crucial role in solving for $X$ because changing the value of $\lambda$ gives a different estimate for $X$.

using $(4),$ is it not possible to auto compute the optimal value of $\lambda$?

If it is possible to auto compute optimal $\lambda,$ why is it in many image de-noising problems we specify the optimal lambda $\lambda$ as an input to the algorithm?

It it is not possible to auto compute $\lambda,$ why is it so?

Thanks for reading my question with patience..


texts: How may samples per symbol in ieee 802.11 g

texts:
I am doing a baseband simulation for IEEE 802.11 g. The standard uses 64 IFFT length ($3.2 \mu s$) and 20 Msamples per second. This result in one sample per symbol. I am wondering what the structure will be if I use 40 Msample per second, i.e will the IFFT size change, or interpolation will be used?

Thanks


texts: What does an image of Fourier Transformation of an image tell us?

texts:
First time studying image processing...
I just don't understand what does fourier transformed image of an image describe?
For example consider given following pictures, The first one is the image, and the second one is the fourier transformation of the image.

Now my question is:
By given the fourier transformation image, what am i suppose to comprehend from that?

I would really appreciate any help here, i just cannot proceed with my studies without understanding this.

Thanks in advance.


texts:
I am trying to understand FFT, DFT through Matlab and I am fairly inexperienced in all of these. I have a time vector t and a corresponding simple sine wave x = 0.25*sin(t). Now I would like to view the 1-sided FFT of this

%MATLAB
%define time and signal
delta_t = 2.6;
t = 0:delta_t:600;
x=0.25*sin(t);
L=length(x);

%define frequency vector
f_low  =1/t(end); %lowest frequency for this range
fs = 1/delta_t; %sampling frequency
f=[1:L]*f_low;

%calculate FFT
FFT = fft(x);

%double sided spectrum
P1=abs(FFT)/L;

%single sided spectrum
P2=P1(1:L/2+1);
P2(2:end-1)=2*P2(2:end-1);
f2=f(1:L/2+1);

%plot
figure()
plot(f2,P2,'LineWidth',1.25);
grid on
xlabel('f')
ylabel('abs(FFT)')

The result of this varies with delta_t (all the while respecting Nyquist theorem). The amplitude of the FFT is different for differnt delta_t and ideally shouldn't it be 0.25, equal to the amplitude of the wave?

Can someone please explain me why does the amplitude vary with delta_t?

Thank you!


texts:
can anybody explain to me what does structure matrix exactly do in image processing? What happens to the pixels?

By structure matrix I mean:

$$S_w[p] = \begin{pmatrix} \sum_{r}w(r)[I_x(p-r)]^2 & \sum_{r}w(r)[I_x(p-r) I_y(p-r)]  \\ \sum_{r}w(r)[I_x(p-r) I_y(p-r)] & \sum_{r}w(r)[I_y(p-r)]^2 \end{pmatrix}$$


texts: Non-separable, higher-order generalized cosine windows in 2D for images

texts:
Higher-order generalized cosine windows define a useful parametric family of windows functions, embedding Hamming, Hann, Blackman, Nutall, etc.:
$$w(n) = \sum_{k = 0}^{K} a_k\; \cos\left( \frac{2 \pi k n}{N} \right)$$
Are there non-separable 2D similar forms, i.e. not a tensor product of the 1D form, usable on the standard rectangular lattice (excluding designs such as in Multidimensional windows over arbitrary lattices and their application to FIR filter design)?


texts: Find the maximum frequency in the signal section?

texts:
I'm recording the received signal after sending a 20Khz signal. I make some rotations with hand before the device and get some variations in the signal.I obtain a very noisy signal and do some pre-processing to it. I tried getting the wavelet response of the signal and this seemded to be good.

I can get the count of the number of times I wave my hand before the device using this. I want to find the change in frequency for every section, the wave of hand happend. I tried applying STFT to the signal.

The response I get is very weak. Is there someway to tweak the signal to obtain the doppler shift ?


texts: Determining sampling rate for signal reconstruction

texts: I want to know how to determine the signal sampling rate required to reconstruct electrical signals whose frequency components are unknown, such as, for example, signals of charge/discharge voltage on a capacitor. I want to design a system capable of scanning this type of signal in real time.


texts:
For discrete time series $Y$, first differentiation ($D_i=Y_i-Y_{i+1}$) and integrator ($S_i=Y_i+Y_{i+1}$) can be defined as two highpass and lowpass LTI digital filters. Where, transfer function for the first difference is:
$$
H(z)=1-z^{-1}
$$
and for the first integrator is:
$$
H(z)=1+z^{-1}
$$
These two transfer function corresponds to $\text{Conv}(Y,[1,-1])$ and $\text{Conv}(Y,[1,1])$ for differentiator and integrator, respectively.

I was wondering what the power transfer function of these two filters should be.

I have to add, in terms of implementation in Matlab, I can get $S_i$ and $D_i$ of a white noise by:

Y = randn(1,100);
D = diff(Y);
S = Y(1:end-1)+Y(2:end);

By filter in Matlab:

S_Filt = filter([1  1],1,Y);
D_Filt = filter([1 -1],1,Y);

Or by convolving conv the corresponding windows:

S_Conv = conv(Y,[1  1]);
D_Conv  = conv(Y,[1 -1]);

Note, all D* and S* variables should be identical (regardless of head and tail of time series, of course).

Thanks in advance.


texts: Inverse DFT in terms of original signals for product of two DFTs

texts:
Given $x[n]$ and $h[n]$ that are zero for $n<0$ and $n>15$, their linear convolution being $y[n]=x[n]*h[n]$ and their discrete time Fourier transforms are $X(e^{j\omega})$ and $H(e^{j\omega})$, then
$Z[k]=X(e^{j\omega})H(e^{j\omega})|_{\omega=\frac{4\pi k}{31}}, k=[0,...,30]$.

How can $z[n]$ be expressed in terms of $y[n]$, where $z[n]$ is 31-point inverse DFT of $Z[k]$.


texts: Incorporating Radial Distortion in Imaging Model

texts:
I'm currently developing a multi-view stereo system but I'm confused as to how all the standard equations in Structure from Motion are affected in the presence of radial and decentering distortion. For example, the standard imaging model is that of

a) World to Camera Coordinates --> translate to camera center and rotate to match orientation of camera (R,t)

b) Perspective Projection --> x = fX/Z and y = fY/Z

c) Conversion to Pixel Coordinates --> scale by no. of pixels/m and add principal point
(b and c can be combined into matrix K)

As a result, we can build a projection matrix P = K[R|t]. Thus a point X in world coordinates can be represented in pixel coordinates as u = PX.

We can also create a fundamental matrix F between two images such that a pair of corresponding point (f on image I1 and f' on image I2) satisfy the epipolar constraint;

x'Fx = 0.

The matrix F depends upon the matrix P. But, my question is that if we incorporate radial and decentering distortion, the matrix P can no longer represent the imaging model. How do we get an epipolar "curve" representation ?

Moreover, for things like triagulation back to 3D coordinates from 2D point.. how do we do it ?


texts:
Find the inverse Z transform:

I have done the solution but my answer does not match with the one given in the textbook. What I may have done wrong?


texts: Segmenting a frequency series: do we need band-pass filtering?

texts:
I am trying to implement what's called the "second spectrum". Basically, you do this:

Take a time series of length $N$.
Divide it into $m$ segment, each of length $N'=N/m$.
For each segment $m$, do a Fourier Transform. The result is the 'first
spectrum',$S_i^{(1)}(f_1), i=1:m$
Divide each spectrum into $n$ octaves, an octave starts at $f_L=f_0 \times 2^p$ and ends at $f_H=f_0 \times 2^{p+1}$, where $p=0,1,2, etc$, and $f_0$ is the lowest frequency in $S_i^{(1)}(f_1)$. For $f_0 = 1$ Hz, The octaves will be like: 1-2 Hz, 2-4 Hz, 4-8 Hz, 8-16 Hz, etc.
For each octave in each spectrum, sum all spectrum values in that octave.
Construct a time-series for each octave. There will be $n$ time series each of length $m$.
Take the Fourier Transform of these time series. This is the second spectrum, $S_n^{(2)}(f_2)$

I have already implemented this algorithm (in c++), but I have a question:

Q: in step 4, can I simply divide the spectrum signal into octaves "just like that"? I mean, without any kind of band-pass filtering to these octaves?

EDIT: what I mean by the "just like that" is as if I am passing the spectrum through a rectangular window (in freq. domain).


texts:
I'm trying to write a correlation function on an embedded device. Because its embedded, I need to make sure that the memory/CPU usage is minimum.

Here's the deal:
- I have a saved waveform sampled at 1 MHz
- I have timing info for the waveform
- I am going to be taking a new waveform at 200 kHz with new timing info
- The timing info for both waveforms is comparable, neither of them start at zero
- I want to determine the delay between one waveform and the next precise to 1 us.

The method I was going to use to do this was:
- Find the maximum of each waveform (there should only be 1 peak in each, and they should match up)
- Set the indexes of each array so that the maximums of the waveforms match up. I'll save the amount of clock cycles that I am shifting in order to do this.
- Then I will correlate this and shift a few samples in each direction and correlate as well. The highest correlation value corresponds to the correct amount of delay, and I will keep this value.

The problem is that my array manipulation skills are limited. I'm struggling when trying to implement this algorithm, especially because one of the waveforms is sampled at 5X the other. Does anybody have any references for this or knows how I can achieve this easily?

The specific problems I'm running into are:
- How do I set the indexes differently so that the maximums match up? I don't want to write a new register because this takes up memory and takes a long time in terms of read/write cycles.
- How do I compensate for the different sampling frequencies? I could fill the 200 kHz waveform with zeros (which won't affect the correlation) but this will take up more memory with no real purpose.
- The timing info isn't precise, in that they don't necessarily line up perfectly. For example, the timing numbers "10000, 10020, 10040, 10060, 10080, and 10100" for the 1MHz signal  doesn't necessarily match up with "10000, 10100" for the 200 kHz signal. So I can't rely on using the timing numbers as an index very easily.

If what I am asking isn't very clear, I can try to embellish a little on requests, but I don't want this to be too long.


texts:
In my systems and signals course I had been asked a question about finding the convolution of two sets.
I was given:
\begin{align}
x[n] &= \{3,2,1\}\\
h[n] &= \{1,-2,3\}\\
\text{Find}\quad y[n] & = x[n]\star h[n]
\end{align}
I have no idea on how to find the convolution of two sets. Any help is appreciated.


texts:
I am new to wavelets and for several times I stumbled over these 2D wavelet transformation pictures, but I do not understand how they are created.

Can anybody try to explain in his own words how this wavelet transformation in 2D works and how or what we can see it on the picture?

I took this image from http://graphicon.ru/oldgr/courses/cg_el01/lect_wav/wav_cg/fig1.gif


texts: DFT of discrete signals, why do we only analyze frequency bins equal to number of input samples?

texts:
If we have a signal $x[n]$ such that we have $N$ samples i.e. $n=0, \ldots, N-1$, then when we analyze the DFT $X[k]$ we only analyze for $k=0,\dots,N-1$ as well.

Why is the range of $k$ tied to the range of $n$?
Why don't we analyze more, or less number of frequencies?


texts:
I have applied transformation(Constant Q Transformation) to my audio time domain signal, this gives me Frequency vs Time representation. Now, the CQT values consists of complex values(84*260), but python makes use of only magnitude values to plot. I saved obtained spectrogram as .png image, to do some image processing work.
So when i read this saved image I get 224*341*3.

After Image processing work, Now I want to reconstruct back my audio time domain signal to check my work.
To reconstruct back to audio domain time signal, first i will read the image, when i read my image I get 224*341*3 values.

I have the phase information of 84*260 resulted from CQT in the beginning. I know we have to make use of phase information to reconstruct the signal. S=Magnitude * Phase.

How would one really reconstruct back the audio signal with this criteria.

Am i doing anything wrong. Has people came across same situations?

I am using Python and librosa library for all the operations.

Edit:

I know librosa icqt has existing bug. Even if I do reconstruction in matlab, how would i tackle the different dimensions. Hope its clear now.

Code:

#load the audio
data_aud,sr_aud  = librosa.core.load(path)

#find the cqt
cqt = librosa.core.cqt(data_aud, sr=44100, hop_length=512, n_bins=84, bins_per_octave=12, tuning=None, norm=1, sparsity=0.01, window='hann', scale=True, pad_mode='reflect')

#take the phase information
mag, phase = librosa.core.magphase(cqt,power = 1)

#save the fig
a_add=librosa.display.specshow(cqt, sr= 44100)
plt.axis('off')
plt.savefig(path,transparent = True, bbox_inches = 'tight', pad_inches = 0)

#image processing
img = cv2.imread(path,1)
gray=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

#multiply the image data with phase????
#reconstruct
y =librosa.core.icqt(imageData, sr=44100, hop_length=512, fmin=None, bins_per_octave=12, tuning=None, filter_scale=1, norm=1, sparsity=0.01, window='hann', scale=True, amin=1e-06)


texts: Is DSP a branch of Computational Science?

texts:
The way I understand it, computational science covers a wide range of fields, all of which use computers and mathematical models to solve a problem involving one of the sciences (physical, life etc.).

But is DSP, which does use maths to manipulate a signal, a form of computational science?


texts: which are saturation and intensity images in this figure?

texts: The attached article talks about hue,saturation and intensity images. Where is each of them in fig6.25 ? I know that fig(b) shows Hue image but where are other two(saturation and intensity images)


texts:
I think this is an easy question but I am confused about the answer. Assume we have correlated noise distributed as $$n \sim \mathcal{N}(0,\sigma^2 \:A),$$ where $A$ is a given square matrix. If we multiply this noise by $A^{-1}$, i.e., $$w = A^{-1} n.$$ Then, the noise covariance matrix in $w$ is given as $$\mathbb{Var}\{w\} = \mathbb{Var}\{A^{-1} n\} = (A^{-1})^2 \mathbb{Var}\{n\} = (A^{-1})^2  A = A^{-1}.$$

Am I missing something?


texts: How to simplify an exponential functions in MATLAB?

texts:
I have many defined functions as below and trying to simplify it with MATLAB simplify() function
$g(t) = 14.9je^{j31°} e^{j2t} - 14.9je^{-j31°} e^{-j2t}$
How can we simplify these kind of functions using MATLAB?
Thank you.


texts: Is there a clever way to implement cascaded moving average filters?

texts:
I'd like to effeciently approximate a gaussian filter's step response using cascaded moving average filters.

I know about recursive moving average, but is there some clever algorithm to cascade them apart from just putting them in series? CIC filters seem related but are mostly used for decimation, I'm interested in pulse shaping.


texts: Intuitive explanation of the Fourier Transform for some of the functions

texts:
Does anyone have a mechanism to understand intuitively  (and automatically) why the Fourier Transform of certain functions have certain shapes (at least for some functions, not necessarily for all)? I know what kind of operator the Fourier Transform is and what it does to a function but somehow I can't see intuitively and automatically why why the Fourier Transform of certain functions have certain shapes. For example, is there a intuitive reason why the Fourier Transform of a pulse (box function) is in the $\mathbb{sinc}$?

Thanks in advance.


texts:
I am calculating spectrogram of a audio file of 36 second using the following code snippet:

window = np.hamming(window_width)
length = signal.size
count = 0
for i in range(0, length-window_width, window_overlap):
    cur_frame = signal[i:i+window_width]
    window_frame = cur_frame*window
    spec = sp.fft(window_frame)
    if(i == 0):
        width = int(length/window_overlap)
        height = int(spec.size)
        specgram = np.empty([height, width])
        energy = np.empty(width)
    energy[count] = getEnergy(spec)
    for j in range(height):
        val = spec[j]
        specgram[j, count] = math.sqrt(math.pow(val.real, 2)+math.pow(val.imag, 2))
    count = count+1
print abs(specgram).max()

I get a ndarray of shape (4000L, 1008L)

but when i calculate spectrogram using Pylab:

Pxx, freqs, bins, im=pylab.specgram(sound_info, Fs=frame_rate)

I get ndarray of dimension: (129L, 125999L).

Actually I need to get the spectrogram which is the dimension of (frequency,nframes)

Thank you in advance.


texts:
Can anybody tell me how to calculate the order analysis / spectrum of a signal? I've been searching in vain. I want to use MATLAB to analyse the signal / sound of a motor.

For more information about the order analysis - here is an example

Thanks in advance for any help!


texts:
Suppose you produce a pulse x in discrete samples of duration dt each (e.g. in Matlab) according to a function.
I was given some notes in which it shows that if I want it to let it have a specific power I have to first normalize it by doing

energy = sum((x.^2).*dt) %first calculate the pulse energy
x = x./sqrt(energy) %normalize its energy to 1

then you calculate the target energy

power = 0 %target power in dBm
power = (10^(power/10))/1000 %target power in W
energy = power*len(x)*dt %target energy of the pulse

and finally you scale the pulse so that it has the target power

x = x .* sqrt(energy) %the pulse at the target power

My question is, how does that really work? Why dividing x by sqrt(E) normalize it to unit energy? And why multiplying the resulting normalized x by the new sqrt(E) brings it to the target power? Why not just divide by E for example? and then multiply by E? Does it have any relation to RMS values? And finally is it something that can be used in any arbitrary pulses, e.g. in sinusoids as well?


texts: Finite and Infinite support in time and frequency domain

texts:
It's known that signals with a finite support in time has an infinite support in frequency (talking about Fourier Transforms).

A common example that is given for the claim is the triangle function with it's Fourier transform $\mathrm{sinc}^2$ function.

I wanted to understand why the following example is not a good one:

$x(t) = \textrm{rectangular}$

$X^f(\omega)=\mathrm{sinc}(\omega)$

We were told in the class that we can't assume that sinc really have infinite support (and that is what I'm not understanding), so the given signal (rectangular) won't be an example for a signal that has a finite support in time and infinite support in frequency.


texts:
In the following lecture:

http://www.ece.mcmaster.ca/~xwu/interp_1.pdf

the model (formula) for solving the linear interpolation problem (1D) given at p.5 is:

$f(x)= a_1x_1 + a_0x_0$

solve for $a1,a0$

the system:

$f(0)=a_i*0 + a_0*1$

and

$f(1)=a_i*1 + a_0*1$

I don't understand where this formula comes from. Even though deriving in all possible ways I could think of, I cannot get the equation of a line $y=mx+b$  ( also written as: $f(x)= mx+b$, where m is the slope and b is the y-intercept) to look something like $f(x)= a_1x_1 + a_0x_0$. Any idea/proposals/discussions etc.?


texts: Establish instantaneous amplitude of signal

texts:
I'm using the YIN algorithm to establish the pitch of my humming / whistling voice.
This is a time domain algorithm and gives me the instantaneous frequency as a function of time ($f(t)$). I now wish to establish the instantaneous amplitude, $a(t)$. I am then plugging these values into the equation:
a(t) * cos(2*pi*integral_of_f(t))

to reproduce the signal based on $f(t)$ and $a(t)$.
At the moment, I am establishing $a(t)$ as follows:
If I establish that $f(t)$ for a given sample is say 20.1 samples in terms of its period, I then create a window centred around this sample t, of length 21. I choose 21 because it is the next whole odd number above the period. I then apply the root mean square equation to this window:
amp = np.sqrt((w * w).sum() / len(w))

My samples are np.float32.
Does this sound like a reasonable approach?
One issue would be if I have a signal with a period drifting between between say 15.01 and 14.99 over time. In that case I would use a window of 17 when the instantaneous period of t is 15.01 and 15 when it is 14.99. Not sure if that is a big issue.


texts: Implementing Welch's method for Power Spectral Density

texts:
I want to implement Welch's method for PSD calculation in MATLAB. I do not want to use built in MATLAB cpsd function, such that I can change the FFT implementation in cpsd(for certain purpose).

I already tried to replicate the method based on Welch paper and explanation from the page: Nonparametric Methods. Below is the code that I already implemented with comments about steps to achieve the result. (it is used for accelerometer real data, so numaxis represent number of axis, $xyz$. Each column represent each axis, the time series data is arranged into row):

function psd = pspectra(data1,data2, window, overlap, nfft, fs)

%checking the size
if(size(data1,2) ~= size(data2,2))
    error('size of column is not equal');
elseif(size(data1,1) ~= size(data2,1))
    error('size of row is not equal');
else
    numaxis = size(data1,2);
end

%calculating window sliding step for iteration
winsize = length(window);
step = winsize - overlap;
iter = 1 + (size(data1,1) - winsize)/step;

%start and end index of first window/segment
istart = 1;
iend = istart + winsize - 1;

%start calculating fft for each window
for i=1:iter
    for j=1:numaxis
        %apply window
        data1(istart:iend,j) = data1(istart:iend,j).*window;
        data2(istart:iend,j) = data2(istart:iend,j).*window;

        %calculate fft
        fft1(:,j,i) = fft(data1(istart:iend,j),nfft);
        fft2(:,j,i) = fft(data2(istart:iend,j),nfft);
    end

    %move to next window segment
    istart = istart + step;
    iend = iend + step;
end

%obtain scale to create modified periodogram
scale = 1/(fs.*sum(window.*window));

%averaging window result and apply the scaling
psd = zeros((nfft/2)+1,size(fft1,2));

for i=1:iter
    psd = psd + fft1(1:(nfft/2)+1,:,i).*conj(fft2(1:(nfft/2)+1,:,i));
end

psd = psd.*scale./iter;

%multiply by 2 except dc and nyquist component (1 and 51)
psd(2:50,:) = psd(2:50,:).*2;

end

I have tried using multiple segment with non overlapping window/segment, and got exactly the same result as cpsd() function from MATLAB. For example, using this call:

pspectra(data1,data1,hamming(50),0,100,100);
pspectra(data1,data1,hamming(25),0,100,100);

However, when I tried overlapping window like this

pspectra(data1,data1,hamming(50),25,100,100);
pspectra(data1,data1,hamming(60),20,100,100);

The resulted power spectra is totally different. I have tried to find if my implementation is wrong (in theory or code), and unfortunately no luck yet. Do I miss something in here?

I really appreciate for any help.


texts:
I know that aliasing occurs when a signal is subsampled. If the sampling rate is lower than twice the max frequency in a signal, aliasing occurs.
How is it in pictures? as far as I know, a sinc-filter is not applicable in photos. so my question is, is aliasing always in pictures?


texts: Implementing Convolution in Frequency Domain?

texts:
Suppose, we have a bitmap image represented as a 2D integer array,
 int [,] image2D; whose FFT is Complex[,] fftImage2D;

Suppose, we have an kernel represented as a 2D integer array,
 int [,] kernel2D; whose FFT is Complex[,] fftKernel2D;

We know that, the convolution (in spatial domain) of image2D and kernel2D would be,

int Rows = image2D.GetLength(0);
int Cols = image2D.GetLength(1);

for(int i=0 ; i<Rows ; i++)
{
    for(int j=0 ; j<Cols ; j++)
    {
        //sweep the kernel2D across image2D
        //...........................
    }
}

We also know that, convolution in frequency domain would be, multiplication between fftImage2D and fftKernel2D.

How can I do this multiplication?

How can I multiply two Complex [,] type 2D arrays of different dimensions? I have understood the theory. My problem is practical implementation. As I described in the question,

Are DFT of the image and DFT of the kernel going to be of different sizes? I guess so. So, how can I multiply them element by element?
In my code, each of the DFTs are represented by 2D Complex numbers. Should, I multiply them according to complex-number's multiplication rule? Probably yes. But, only when their dimensions are same. Right?


texts: Transition table / matrix for convolutional encoder K=3 understandings

texts:
Im trying to understand convolutional encoder $K=3$ , and I almost understand it but still confused about what's called transition matrix which in other words it means

what input value ($0$ or $1$) would produce the next state, given the current state.

So if we have $K=3$ which means $3$ shift registers, there's $2^{K-1}$ possibilities which means there's $4$ states as Im showing down: (it's in other words state machine)

And I fully understand it, and Im stuck on how do I built the correspond transition matrix for this state machine?
I've seen on the internet something weird that the transition table for convolutional encoder $K=3$ is:
transitionMatrix = [0 1 1 55; 0 0 1 1; 55 0 55 1; 55 0 55 1]

it was written that value $55$ is value of STATE HISTORY .
I didn't understand this matrix and how this matrix describe the state machine that I attached the photo above for convolutional encoder $K=3$
Could anyone please help me out for this?
I already did the state machine and I understand it, I didn't understand how do I build the transition matrix. I'm posting here in order to understand it and not just copying it or reading it from a book without knowing how do I implement it.


texts: Getting Frequencies Corresponding to Peaks In FFT Plot - MATLAB

texts:
So I'm currently working on a project where I need to find the frequencies for specific peaks in my FFT plot of an audio file. However, I'm not sure how to actually get the frequencies that correspond to the peaks I'm interested in.

This is my code: https://pastebin.com/zZW8B5Th

And this is a picture of my FFT plot: https://imgur.com/0ajtful

Essentially, I'm trying to construct a Fourier Series from the signal. Although I know how to find the coefficients of the sine/cosine terms, I'm just a little confused on how to find the frequency.

Thanks everyone!


texts: How can I classifying this License Plate images? Matlab , LPR ,

texts:
I want to perform automatic character recognition on licence plates. One of my preprocessing step is binarization. However I have four different class of images and I have to use a different automatic thresolding algorithm on each class to get good binarization results.

I would like to choose the right algorithm automaticaly, my idea is to first perform classification on the licence plates images. In the following picture, each row represent a class.

picture:

http://www.mathworks.com/matlabcentral/answers/uploaded_files/31813/class.jpg

I'm working with MATLAB, how can I do that?? what is the best way?
Any ideas are appreciated.


texts:
This may be a very straightforward question but I've been trying to understand this for a long time and can't wrap my head around it.

How would one go about discretizing the following continuous time signal:

signal_continuous = 1/6 * sin(880*pi*t) + 1/15*sin(1760*pi*t);

As far as I understand, if I wanted to sample the above signal with a sampling rate of 8Khz over two seconds. I would do something like this (in Matlab speak):

t = [0:1/8000:2-1/8000];
signal_discrete = signal_continuous./t;

Is this the whole picture or have I completely missed the point here?

My rationale behind doing the above is that since the continuous time signal is, well, continuous and if I want to sample it at 8000 cycles per second over a two second period. I would just divide the signal by that (t in this case). However, what happens to the rest of the signal? Or rather, does this not distort the whole signal (if the original was longer than 2 seconds, for example)?

EDIT:

I was complicating things for no reason. In Matlab I would just need to do the following:

t = [0:1/8000:2-1/8000];
discrete_signal = 1/6 * sin(880*pi*t) + 1/15*sin(1760*pi*t);

This would create a two second discrete-time signal out of the original continuous signal.

The theory provided in the answer, I think, is beneficial so I'm leaving the question up.


texts: harmonic waves as integer multiple in spectrum

texts:
i have a motor that is rotating with a certain frequency. If i check the frequency spectrum it contains a peak on 150 hz. Also i have peaks at 300,450,600 ... i guess that those peaks are harmonics.
But i dont really understand why harmonics appear.

I think the reason is:
a periodic signal can be described by sine, cosine waves (except the if periodic signal itself is a sine or cosine wave) So because a rotating motor is periodic, the frequency contains multiple cosine or sine waves that together build the signal. And the harmonics are the parts of the signal?

but I still don't understand why they occur if my machine is rotating at 150 hz. How are the other frequencies build?


texts: How to analyze reverb in existing sound recording?

texts:
I'm wondering if there is a technique by which to analyze an existing sound recording and extract some details of the reverb and echo that are present (either due to natural environment, or added during production).

Are there software packages or libraries that implement that sort of analysis? (Is there even a chance of recovering the audio minus the reverb?)


texts: Shifting a fourier spectrum by subpixel amount in python

texts:
I am working on a Fourier Ptychography problem. My research problem requires me to shift the Fourier spectrum of an image by a floating-point value. For a real-valued image, we can simply use cv2.warpAffine to shift the image by floating-point values (aka subpixel shifting).

Taking fourier transform of an image produces a complex matrix. The problem is, cv2.warpAffine does not support complex matrices, and so I cannot use it on them. I tried searching for alternatives, but none of them seem to work. I came across numpy.roll, but the problem is, it does not support subpixel shifting. Rounding off the shifting values translates to loss of information in my case. Is there a solution in python, that allows for subpixel shifting on complex matrices?

Thanks.

EDIT: Based on Marcus' answer, I did some digging and implemented a nifty little script for subpixel shifting in python based on a Matlab script for the same.

Here's the link to the script. Hope it helps!


texts: The Hilbert Transform of an Amplitude Modulated signal returns the envelope of the signal. What does the Hilbert transform of a Frequency Modulated signal return? How can I use the Hilbert Transform to get the sidebands of a Frequency Modulated signal?


texts: Is document image binarization a closed research field

texts: A lot of work has been done in the field of document image binarization. Is it a closed research field or still an interesting challenge area? I am new in the field, and just curious.


texts: Detecting periodic components in discrete-time signals

texts:
I am looking for real data examples where it is of interest to detect the presence of a periodic component of some discrete-time signal when the period of the periodic component is not known.
Suppose that $\{X_t:t\in\mathbb Z\}$ is a discrete-time signal (univariate or possibly multivariate) and consider the situation where this signal can be modeled as $$X_t=s_t+\varepsilon_t$$ for each $t\in\mathbb Z$, where $\{s_t:t\in\mathbb Z\}$ is a deterministic periodic component such that $s_t=s_{t+p}$ for all $t\in\mathbb Z$ and $\sum_{t=1}^ps_t=0$ with some period $p\ge2$ while $\{\varepsilon_t:t\in\mathbb Z\}$ is some stationary random noise sequence. We observe $X_1,\ldots,X_n$ and we want to detect the presence of the deterministic periodic component $\{s_t:t\in\mathbb Z\}$, i.e. we want to know whether $X_t=\varepsilon_t$ for $t\in\mathbb Z$ or $X_t=s_t+\varepsilon_t$ for $t\in\mathbb Z$ with some period $p\ge2$.
In many situations, we expect the presence of certain periodic components and hence it is not reasonable to assume that $p$ is unknown. For example, if we record temperature, then it is clear that we should expect to see daily and yearly periodic components. It does not make a lot of sense to check, say, for a three day periodic component. In many situations, some theory suggests that there should be some periodic components present and then it makes sense to try to detect these particular periodic components instead of assuming that the period is unknown.
However, I suspect that there are examples (maybe in engineering or astronomy) where it is not clear whether the periodic component should be present and what period of the periodic component we should expect to detect.

I am looking for real data examples (preferably multivariate but univariate examples are also of interest) where it is not clear whether the periodic component should be present and where it is not clear what period we should expect to detect. References to textbooks, papers or websites are very welcome.

I have no background in signal processing and my question or phrasing might seem strange. Please let me know if I should clarify some parts of my question.
Any help is much appreciated!


texts: Filter signal for beamforming and phase response

texts:
I was wondering this for quite some time but couldnt find a satisfied answer. Seems like not a lot worry about phase response in general. But maybe I'm thinking about it in a wrong way...

I want to apply a high pass filter to my microphone pressure signals to cut away the low frequency garbage before beamforming. So I was wondering about using filters (low, high, band) on your time signals before beamforming.

How do you guys correct for the phase response? And is it actually needed? If not, is only applying (e.g.) a butterworth filter adequate enough to perform correct beamforming afterwards (outside the frequency for which it is filtered)? Or is some correction in phase needed?


texts: Accurate measurement of object's position on a bayer-filtered and grayscaled images

texts:
This is an experience-related question about cameras and image processing.

I have a laboratory setup. Under some simplifications it looks like a red ball rolling on a white table (the real setup is shown on the video). A computer vision system should measure X and Y positions of the ball.

For now I use a color camera (with Bayered output) and I find the ball by applying a simple color filter (InRange for R, G, and B components) to the frames. It works, but I would like to increase precision of the measurements.

Question:
Would accuracy of measurements be better in a case I used a grayscaled camera with the same parameters as bayer-camera?

I guess it should be so in a perfect case, i.e. when the ball is definitely highlighted on the frames. But in reality brightness filter usually works worse than color filter. On the other side I can change the shooting conditions, because this is laboratory setup. I can change the color of the objects, I can change lighting.


texts: Filter, window function and audio signal: how to link them to each other?

texts:
I have a question about the implementation of digital filters in CoreAudio. I'm in big trouble because it is a few weeks I'm trying to understand how to implement them. The basic idea is this: while I talk to the iPhone's microphone, my voice is filtered by a low-pass filter or high-pass or bandpass.

Reading the book "Digital Signal Processing" and following the algorithms mentioned in the book "Digital Sound & Music" I came to write the following code (high pass filter):

float cutoff = 5000;
float rate = 44100;
int N = 1024;
float f_c = cutoff/rate;
float ω_c = 2*M_PI*f_c;
int middle = N/2;
int i = -N/2;
float fltr[N];

do {

    if (i == 0)  {
        fltr[middle] = 1 - 2*f_c;
    }
    else {
        fltr[i+middle] = -sin(ω_c*i)/(M_PI*i);
    }
    i++;

} while(i != N/2);

After writing the filter, in the book, you are instructed to run a windowing function. Here I have the problem: I do not understand what to do. Reading on Apple's official documentation, I find the function of Blackman.

int N = 1024;

float dataOut[1024];

for (int n=0; n < 1024; ++n)
{
    dataOut[n] = 0.42 - (0.5 * cos(  2 * M_PI * n / N ) ) + (0.08 * cos( 4 * M_PI * n / N) );
}

But once inserted, what should I do? How can I connect the filter carrier, the Blackman function and the audio signal? I just can not understand it.

Thank you for your attention.

UPDATE 1

I answer to users bristow robert-johnson and Matt L. (I write here because the comment was too long to be invited)

Hi, thank you for answering. There is no filtering at this time because I do not understand how to connect the filter (?) With Blackman window.

What I understand from the chapters devoted to the book “Digital Signal Processing” filters: Create a filter, a sort of kernel filter that is a kind of "mask" to be applied to the signal itself. In this state the filter is said "Impulse Response" or also is said IIR. Now I apply the mask to signal by a function convolution becoming a FIR filter. What I have written is this correct?

As for calculating coefficients, as I understand it, they should be the coefficients Nyquist angular frequency. Is right? If yes, should be variables of "f_c" and "ω_c" (this is written in the commentary of the algorithms in the book).


texts: FFT of some signal $x$ has two parts: $|x|$, which is the magnitude spectrum and $arg(x)$, which is the phase spectrum. I want each frequency to have the same phase. Is there a way to accomplish this?


texts: Why is CDMA so complicated, wouldn't this simpler solution work?

texts:
I am now reading about Code Division Multiple Access method (CDMA), which is used for multiple hosts to send data at the same time.

The main idea behind it is to use orthogonal vector basis to encode multiple bit streams into one, send the joined stream across the link, then decode the stream differently for every separate message.

What confuses me is that in digital logic negative numbers are introduced, as well as 0 indicating no signal.

When I first heard that such encoding exists, my first idea of how it should work was different from the original and it looks way simpler.

Could someone check and add comments to my method? WHY is it not used instead of ORIGINAL CDMA?

Outline of my solution

I take n orthogonal vectors (where xor is the addition operation in my vector space), each vector corresponding to one chip code. I then encode every bit as a chip code if it is 1, as 0 otherwise. Finally, I xor all bit sequences together and send the joined result. At the other end, recipient can decompose the bit stream into unique components and read the message of every host.

Solution in more details

Assumptions

There are n hosts (numbered from 0 up to n - 1), and all of them want to send b bits at the same time using the same link.
The receiver gets a joined bit stream, but is interested in reading the message of only one host (WLOG the message sent from the 0-th host).

Problem: create a protocol that would satisfy the assumptions.
My solution

For i-th host, assign a chip code to him of the form 2^i (power not xor here) written in binary. Call this chip code c(i). For example:

Chip code 1000..00 would be given to host 0
Chip code 0100..00 would be given to host 1
...
Chip code 0000..01 would be given to host n-1.

Every single bit will be encoded using n bits.

if (host i wants to send 0) {
    encode it as 0000...00  // n zeroes
} else {
    encode it as c(i)
}

For host i wishing to send j-th bit, call the encoded bit as e(i, j). For example:

if (host i wants to send a bit stream, such that j-th bit is 1) {
    e(i, j) = c(i)
} else {
    e(i, j) = 0;
}

For every j, compute

d(j) = e(0, j) ^ e(1, j) ^ ... ^ e(n - 1, j)  // xor operation here

Concatenate d(0), d(1), ... d(n - 1) and send these n * b bits through the link. This will be the joined stream that the recipient at the other end should receive.
The recipient wishing to read the data sent by host 0 can look at (0-th, n-th, 2n-th, ..., b(n-1)-th) bits to see what the first host has sent.
The solution is simpler than original, because it uses only xor operation for encoding, and digits are only binary (i.e. no -2, +3 when doing vector addition). QED.

Does this solution work? If yes, why is it not used? The best reason I could guess is that we want to distinguish between hosts sending something and nothing. Are there other reasons?

Aside question: what are "nearly orthogonal" vectors in the context of CDMA?


texts:
Let's say the task is to determine element position on image.
First very important thing is correct detection of object then some algorithms of calculating position are used (for examble blob analysis). Everything depends on multiple things (detection correctness, used algorithms etc.)

Lets assume we have callibrated image and know error given by callibration.
What are the methods to calculate reliably precision of computer (and machine) vision algorithms?
Can it be done analiticaly or only by experiments and tests?

The question addres cases when we detect element position and also other computer vision problems.

I want to get references to problems which are related to computer/machine vision especially element position detection and present some correctness computations either analyticaly or experimental approach to show this precission.

Also suggestions how to improve this question are welcomed.


texts: What is the advantage of using a Stepped Frequency waveform over FMCW?

texts: What is the advantage of using a Stepped Frequency Waveform instead of a traditional linear FMCW ?


texts: filterbank: understand the different responses at the center frequency of each filter

texts:
I'm computing filterbank by applying 26 triangular filters on a Mel-scale to the power spectrum of an audio frame to extract frequency bands, and I found that some references use filterbank with equal responses at the center frequency of each filter, see the following figure:

A formula for calculating these is as follows:

Then, they take the log of the resulting energies to compute the log filter bank energies.

But, in some other references, they use filters with decrease responses at the center frequency of each filter, i.e. something like this:

I want to understand the difference between these two methods, and which one is more accurate for usage in a speech recognition application based on deep learning?

And what are the formula and different steps to calculate the second one?

References:

http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/

http://www.cs.cmu.edu/afs/cs/user/bhiksha/WWW/courses/yahoo2009/01-02.featurecomputation.ppt


texts: Simulink Floating-Point and Fixed-Point Modeling

texts: I have not used Simulink extensively. I would like to know how to determine if a model is fixed-point or floating-point. Also, how to convert from fixed-point [back] to floating-point. The posts on the internet mostly discuss floating-point to fixed-point conversion. I need fixed-point to floating-point to learn and be able to use both floating-point and fixed-point models.


texts: How did Scientists Capture First Images of Molecules Before and After Reaction?

texts:
"Using a state-of-the-art atomic force microscope, the scientists have taken the first atom-by-atom pictures, including images of the chemical bonds between atoms, clearly depicting how a molecule's structure changed during a reaction. Until now, scientists have only been able to infer this type of information from spectroscopic analysis."

http://www.sciencedaily.com/releases/2013/05/130530142007.htm

Are there any papers how this thing actually works?


texts: Q matrix and updating times in a Kalman filter

texts:
The context of the problem is that I have several robots located remotely which give their position (x,y coordinates) every x seconds and send it to a centralized remote server. The value of the seconds, x, depends on the update frequency of the robot, so it can be 60, 70, 80 or whatever seconds.
In this server, I take these measurements (x, y coordinates of each robot) and generate a Kalman filter to predict next steps, since the measurements are given each "too much" time and I need these predicted positions.
The Kalman filter has the equations:

x[k] = x[k-1] + v[k-1]*dt
v[k] = v[k-1]

The robots are linked by pairs, and the filter of each robot (inside a pair) is updated when a measurement arrives from any of the two robots.
Here comes the problem. The frequency of update in the filters depends on the time of the robots' measurements. For instance:

Two robots giving measurements every 10s and every 25s. The filters will be initialized at t=0s and then updated at t=10s, t=20, t=25s, t=30s, etc.

Two robots giving measurements every 3s and every 5s. The filters will be initialized at t=0s and then updated at t=3s, t=5, t=6s, t=9s, t=10s, etc.

Because of this different updating times, the filters of some pairs of robots work very well, but others do not.
The matrixes Q and R stays constant at any moment. For the Q matrix, the values where given acording to the error produced by the maximum aceleration that can be reached (the problem should be here).
How can I implement a filter (a class) that works for all of the pairs independently of the updating filter times and updated measurement times? Which is the dependence between the update time of the filter and the time when a measurement arrives? I mean, how does it affect quantitatively?
Should I modify the matrixes Q, R and/or P according to the time diference?
It is an open question since I am a bit lost. Any help would be highly appreciated.


texts: Help generating coefficients for basic FIR LPF

texts:
I'm trying to generate coefficients for a FIR low pass filter to be used in an FPGA.  I'm using python (scipy.signal) to attempt to do this, but am having trouble getting coefficients in a usable form.

This is the function I'm using to generate coefficients:

b = signal.firwin(21, 0.01, window='blackman');

The coefficients it returns are:

[ -1.62779573e-18   1.08170091e-03   4.74493669e-03   1.19927580e-02
   2.37996986e-02   4.03774280e-02   6.06306778e-02   8.20599398e-02
   1.01201426e-01   1.14487974e-01   1.19246920e-01   1.14487974e-01
   1.01201426e-01   8.20599398e-02   6.06306778e-02   4.03774280e-02
   2.37996986e-02   1.19927580e-02   4.74493669e-03   1.08170091e-03
  -1.62779573e-18]

These are unusable to me when it comes to FPGA implementation.  How do you scale these or is there a better technique that produces either integer coefficients or coefficients that could be scaled to fixed point in hardware?


texts: Can the Entropy Be Used as a Measure of the Contrast of an Image?

texts:
I was looking for value that can measure a contrast of an image. And I found a couple of answers like this one https://stackoverflow.com/questions/13397394/how-to-measure-contrast-in-opencv-visual-c/13398126.

But I can't find any explanations WHY is that. No articles or links.
So I wrote this code:

def shannon_entropy(img):

    histogram = img.histogram()
    histogram_length = sum(histogram)

    samples_probability = [float(h) / histogram_length for h in histogram]

    return -sum([p * math.log(p, 2) for p in samples_probability if p != 0])

And found the Shannon entropy for four images:

 - original image "image_1". Entropy is 7.82426730799067

 - more contrast image "contrast". Entropy is 5.682661230050924

 - much more contrast image "contrast2". Entropy is 2.1507540571698587

 - less contrast image "contrastLow". Entropy is 6.565885023104716

I thought, that for the last image entropy will be the highest, but no. And it completely baffled me.

Can somebody explain:

Why entropy is a measure of contrast?
What's wrong with my images or calculations?


texts: How to Pick a Journal for Publishing a Research Work?

texts: I got my MS in EE way back in 1994 and have not presented at a conference, or published in a journal before, but I had and still have an intense interest in the applied math that I learned. I have been working my own research project for a few years, and I want to publish a paper on my work. What I did is develop an algorithm similar to a discrete Kalman fiter, and it gives better estimates under certain conditions. I am looking for advice on what journal would a good choice. It would take to long for me to learn LaTeX (is that what it's called), and I don't have the education needed to develop the sophisticated math that I see in many journals.


texts: Features of a non-periodic signal for comparison

texts:
I am looking for comparing a signal (256 / 512 point discrete | non-periodic ) (obtained after extracting some features from the captured image) with a database of around 5000 signals (256 /512 point discrete | non-periodic)(constructed with the same features as the one in the image). The features are angles and distances from an image captured. The signal is constructed by taking angular value on the x-axis and the corresponding radial distance on the y-axis (amplitude). See the two images.

This image has 15 non-zero points and can be considered as one of the 5000 entries in the database.

Second image - signal acquired:

This image is that of the signal acquired and as can be seen it has 12 non-zero points in the signal. This is actually the corresponding signal of the first image but with some missing points and some variations in the amplitude.

The above is just a example of one of the entries of the database containing 15 non-zeroes and the acquired signal containing 12 non-zeroes. There can be many entries in the database containing the same number of non-zeroes, but the positioning of the non-zeroes will, of course, be different. The number of non-zeroes in a particular entry of the database varies from 5 to 50.

Now consider this as a signal to be compared.

There are two bottlenecks - Firstly, the acquired signal contains missing points from its corresponding signal in the database. The missing points can be very high sometimes. This can be viewed as the correlation between two signals can drop from 1 to 0.5 sometimes. Still I can apply correlation and look for the correct signal match.

But, here comes the second bottleneck and that is I do not want to correlate the acquired signal with each and every entry of the database (because it takes a huge amount of time). Instead, I want to shortlist some entries first based on some feature that can be calculated from the signal. This feature can be mean, variance, standard deviation, 3 dB bandwidth, etc. However, this feature must be something that does not vary much in the case that some points are missed as I stated earlier. I also thought of applying FFT on the signal and using some feature such as mean of the absolute value, phase information, etc.

Eventually, I have tried mean, variance, std, power bandwidth, mean of the absolute value of FFT. However, the change is quite drastic and I cannot shortlist to a small number (by small number I mean around 1000). At the same time, as I mentioned earlier, the correlation also drops quite drastically in some cases thus failing to find a match. Do you have any suggestions on how to go about for this problem of comparison when the acquired signal has missing points from the corresponding database signal? Any idea is appreciated! Thanks!


texts: Reducing the number of data points of an FFT

texts:
I'm in an odd situation where an FFT is giving me too much data to process. I'm performing an FFT on audio sampled at 500 Hz, over 30 seconds. The result is 7500 bins of useful information, however I don't need this sort of resolution. I still need to see frequencies up to 250 Hz so I can't reduce the sampling rate. I also need information from the entire sample, so I can't use less than 30 seconds of data either.

What can I do to reduce the number of points, but keep the magnitude and phase information relatively intact? I would prefer not to throw away points, is it possible to take an average of x points?


texts:
I try to remove from an image all stripes with slope between 110° and 120° (Fig. 1a) and my first idea was to filter my image in the frequency domain.

As far as I remembered, the spectrum of those stripes should occupy a pie segment between 20-30° (i.e. the original slope - 90°), but to be sure I prepared an artificial black-and-white image with similar pattern (Fig. 2a). It seemed my speculations were correct (Fig. 2b).

Now, I tried to create a mask and simply "cut off" part of the spectrum to supress high frequencies associated with 10-120° stripes.

Only after several trials I came across a mask (Fig. 3a) which satisfies my needs to some degree (Fig. 4). However, I can't figure the relation between this empirically-obtained mask and the spectrum.

Could someone give me some advice where are flaws in my code or in my reasoning?

Here's the MATLAB code I've used for filtering:

f = imread('img/crystal-grain_sample.jpg');
figure, imshow(f); title('f');

m = max(size(f));
P = 2^nextpow2(2*m);
PQ = [P, P];

Fp = fft2(f, PQ(1), PQ(2));

Fc = fftshift(Fp);
Fc_abs = abs(Fc);
Fs_log = log(1+Fc_abs);

%figure; imshow(Fc_abs, []); title('Fc (abs)');
%figure; imshow(Fs_log, []); title('Fs (log)');

Hp = im2double(imread('spec-mask.png'));
Gp = Hp .* Fp;
gp = real(ifft2(Gp));
gpc1 = gp(1:size(f,1), 1:size(f,2));
figure; imshow(gpc1, []); title('g');

And here's the original image:


texts: Playing back an audio file at an arbitrary rate with upsampling and downsampling

texts:
Background: My overall goal is to create an audio file player which can playback an audio file at an arbitrary rate. The rate can change continuously as the file plays. This means upsampling to slowdown the playback when the rate < 1, and downsampling to speedup the playback when the rate > 1. (I'm not concerned with maintaining pitch.)

Let's focus on the slowdown (upsample) part to make the question simpler.

I read that to upsample by a rational fraction L/M, I should first

Upsample by L
Apply an interpolation (lowpass) filter
Downsample by M

However, this doesn't work well for an online, live context, because L or M can be very large depending on the desired playback rate, requiring large buffers for upsampling/downsampling. I want to process the audio file in small chunks as it plays, so my buffers must be kept to a reasonable size.

I decided to try the following algorithm instead, for upsampling by an arbitrary floating point rate:

// pseudo-code upsample algorithm
sloppy_upsample(int source[], int destination[], float rate, int frames):
  for each frame i:
    destination[i] = source[round(i * rate)];

This is a sloppy "zero order hold" algorithm that probably doesn't get the sample rate exactly correct, but it actually sounds surprisingly good for my purposes. I also like this algorithm for its simplicity. But, of course, there's audible aliasing.

I know that I can apply a lowpass filter before upsampling to prevent aliasing, but I'm not sure what frequency to set the cutoff at. In similar situations such as this one, the cutoff could be set to less than L * 0.5.

Is it possible to select a good cutoff for my sloppy_upsample function above that prevents aliasing?
Is there another approach for upsampling/downsampling that works well for live playback?


texts: What effect does a delay in the time domain have in the frequency domain?

texts:
If I have a signal that is time limited, say a sinusoid that only lasts for $T$ seconds, and I take the FFT of that signal, I see the frequency response. In the example this would be a spike at the sinusoid's main frequency.

Now, say I take that same time signal and delay it by some time constant and then take the FFT, how do things change? Is the FFT able to represent that time delay?

I recognize that a time delay represents a $\exp(-j\omega t)$ change in the frequency domain, but I'm having a hard time determining what that actually means.

Practically speaking, is the frequency domain an appropriate place to determine the time delay between various signals?


texts: low frequency issue with noise spectrum (of derivative)

texts:
I have hard time in understanding the results I get when I plot the spectrum of the time-derivative of a noise signal:

$|DFT[\frac{dx(t)}{dt}](\omega)|^2$

as a function of frequency $\omega$ (see image below), where $DFT$ is the discrete Fourier Transform, $x(t)$ is a discrete time serie of random samples extracted from a gaussian distribution.

In python:

import numpy as np
import matplotlib.pyplot as plt
x = np.random.randn(100000)
F = np.abs(fft.fft(np.diff(x)))**2
plt.loglog(F[:50000])

The F.transform of a derivative is $i\omega$ times the F.transform. So I'd expect the result to be  $\omega^2 |DFT[x]|^2$.

This happens, but what are the features at low frequency? If I run the script several times, a plateau shows up at low frequency, but not always with the same amplitude and not always at the same low frequencies (see in the fig. the differences between various runs of the same script). Probably something strange happens also at high frequency? I heard "leakage" can explain things like this, but I'n not an expert. Can you help me to understand, and correct it?


texts: Where is the Nyquist frequency in an FFT spectrum?

texts:
I've been working with even length audio for a while, so I've had no problem finding the Nyquist frequency in the FFT spectrum, as it's shape is:
$$(\text{dc_component}, f_1, f_2, ..., f_{N/2-1}, f_\text{NYQUIS}, -f_{N/2-1}, -f_{N/2-1}, ..., -f_1)$$
However, an odd number of elements doesn't fit this. What's the layout for odd FFT spectrums?


texts: How can I calculate what a recording would sound like if it was played through my speakers and mic?

texts:
I am trying to isolate the sound of someone speaking while music plays back in a noisy environment.

I have multiple speakers and a microphone in my system. I am able to

play back just music and record it, then
play back music while speaking and record it.

After that, I can look at the difference in the signals and get the meaningful information I need.

In order for this to work in real-time, I would need to calculate what that first recording would sound like. If I have the impulse+magnitude response or some other data about the system, is it possible to predict what the microphone will detect?

EDIT:I suspect that some kind of transform could be derived from the impulse response that will give a decent prediction of what the mic will hear. To be clear, this doesn't need to be incredibly accurate. I just want a general idea of what the mic will pick up.


texts: What is going wrong with my IDFT program?

texts:
I learn IDFT and have created two programs.

Input x(n). It calculates input x(n) values from 2 sinusoids and outputs them.
Output x(n). It performs IDFT on same signal and outputs results.

Results should agree. Unfortunately while first program's output is okay, second's is not. Can you please find the error and say me?

Input x(n) program

// Input x(n) program.
// Grid is created by variables lowTestFreq and highTestFreq
// both of double type, and integer type M.

#include <iostream>
#include "math.h"

int main()
{
//Unchanged variables.
float pi = 3.141592653589793; //pure number.
int n = 0; // sample.
float input_xn = 0;

//Sampling frequency and samples of segment.
int samplFreq = 100; // sample/second.
int N = 40; //sample.

// 1st existing sinusoidal.
float ampl_1 = 70000; //pure number.
float existFreq_1 = 36; // Hz.
float phase_1 = 0.67 * pi; // rad.
float unitAngle_1 = (existFreq_1 / samplFreq) * (2 * pi); // rad/sample.

// 2nd existing sinusoidal.
float ampl_2 = 600000; // pure number.
float existFreq_2 = 25; // Hz.
float phase_2 = -0.37 * pi; // rad.
float unitAngle_2 = (existFreq_2 / samplFreq) * (2 * pi); // rad/sample.

for (n = 1; n <= N; ++n)
{
input_xn =  ampl_1 * sin(n * unitAngle_1 + phase_1)
      + ampl_2 * sin(n * unitAngle_2 + phase_2);
std::cout << "input-x(n)" << input_xn << std::endl;
}
return 0;
}

Output x(n) (IDFT) program

// Grid is created by variables lowTestFreq and highTestFreq
// both of double type, and integer type M.

#include <iostream>
#include "math.h"

int main()
{
//Unchanged variables.
float pi = 3.141592653589793; //pure number.
int n = 0; // sample.
float input_xn = 0; //pure number.
float Re = 0; //pure number.
float Im = 0; //pure number.
float xn_cos = 0; //pure number.
float xn_sin = 0; //pure number.
float co = 0;
float si = 0;
float output_xn_Re = 0;
float sumReal = 0;

//Sampling frequency and samples of segment.
int samplFreq = 100; // sample/second.
int N = 40; //sample.

// 1st existing sinusoidal.
float ampl_1 = 70000; //pure number.
float existFreq_1 = 36; // Hz.
float phase_1 = 0.67 * pi; // rad.
float unitAngle_1 = (existFreq_1 / samplFreq) * (2 * pi); // rad/sample.

// 2nd existing sinusoidal.
float ampl_2 = 600000; // pure number.
float existFreq_2 = 25; // Hz.
float phase_2 = -0.37 * pi; // rad.
float unitAngle_2 = (existFreq_2 / samplFreq) * (2 * pi); // rad/sample.

// Testing grid.
float testFreq = 0; // Hz.
float testUnitAngle = 0; // rad/sample.
float testAngle = 0; // rad.
float lowTestFreq = 0; //Hz.
float highTestFreq = samplFreq; //Hz.
int M = N; // pure number.
int m = 0; // pure number.
float testFreqStep = (highTestFreq - lowTestFreq) / M; // Hz.

for (m = 0; m <= M-1; ++m)
{
testFreq = lowTestFreq + m * testFreqStep; // Hz.
testUnitAngle = (testFreq / samplFreq) * (2 * pi); // rad/sample.

Re = 0; Im = 0;

for (n = 1; n <= N; ++n)
{
input_xn =  ampl_1 * sin(n * unitAngle_1 + phase_1)
      + ampl_2 * sin(n * unitAngle_2 + phase_2);

testAngle = n * testUnitAngle;
xn_cos = input_xn *  cos(testAngle);
xn_sin = input_xn * -sin(testAngle);

Re += xn_cos;
Im += xn_sin;
}
sumReal = 0;
for (n = 0; n <= N-1; ++n)
{
co = cos(2 * pi * m * n / N);
si = sin(2 * pi * m * n / N);
output_xn_Re = Re * co - Im * si;
sumReal += output_xn_Re;
}
std::cout << "output-x(n): " << sumReal / N << std::endl;
}
return 0;
}


texts: Analyze stability of a closed-loop system with Bode

texts:
I have an open-loop system function $L(s)$ and its Bode plot is

As MATLAB says, it is stable if we close the loop with unitary feedback. I thought that, seeing the Bode plots one could tell if the closed-loop system would be stable if the $0\textrm{ dB}$ crossing occured at a lower frequency than the $-180°$ crossing. I made the Blode plots for $0.01L(s)$ and got the following:

Now the closed-loop system would be stable too, but this time the $0\textrm{ dB}$ crossing occurs at a lower frequency than the $-180°$ crossing. Nevertheless, in both cases the closed-loop system turns out to be stable.

Then I made the Bode plots for $0.1L(s)$ and got this:

And now the closed-loop system is unstable.

So, my question is: how can one know, with just looking at the Bode plots, if the closed-loop system is going to be stable or not?


texts:
I cant figure out why our output image is 257 * 32 , 32 I know why , but 257 I cant understand ? this is the exact word of this paper wich applied STFT on 3 channel of eeg signal :

Short time Fourier transform (STFT) was applied on the
time series for each 2 s long trial. In case of 250 Hz signal this
is corresponding to 500 samples. STFT was performed with
window size equal to 64 and time lapses equal to 14. Starting
from sample 1 toward sample 500, STFT is computed for 32
windows over 498 samples and the last 2 samples, remaining
at the end, are simply ignored. This leads to a 257   ×   32 image
where 257 and 32 are the number of samples along the frequency and time axes respectively.

I really really appreciated if anyone could explain it to me


texts: given an image and its fourier transform, find another image fft using fourier properties

texts:
I managed to solve few of those questions , but I'm not manage to answer 3 of them correctly,

the question is: given an image and its fourier transform, find another image fourier transform ,using fourier properties to explain why it would look like that.


texts: Example of non-equivalence of the two PSD definitions

texts:
According to the book Introduction to Spectral Analysis by P. Stoica and R. Moses, the power spectral density (PSD) $P(\omega)$ can either be defined as the discrete-time Fourier transform (DTFT) of the covariance sequence $r(k)$, i.e.,
\begin{align}
P(\omega)=\sum_{k=-\infty}^{\infty}r(k)e^{-j\omega k},
\end{align}
or alternatively as
\begin{align}
    P(\omega)=\lim_{N\to \infty}E\Bigg\{\frac{1}{N}\Bigg|\sum_{n=0}^{N-1}x(n)e^{-j\omega n}\Bigg|^2\Bigg\},
\end{align}
and these two definitions are equivalent under the assumption that
\begin{align}
    \lim_{N\to \infty}\frac{1}{N}\sum_{k=-(N-1)}^{N-1}|k||r(k)|=0.
\end{align}
What is an example of a signal for which these two definitions yield different results, i.e., for which the above assumption does not hold?


texts: EWMA based differential that is strictly positive on strictly increasing irregular time series?

texts:
Let's say I have a time series $ \{t_i, x_i\} $ where the $x_i$ are strictly increasing and the $t_i$ are not regularly spaced. I can define a differential operator as

$$ \Delta_i = \frac{x_i-x_{i-1}}{t_i-t_{i-1}} $$

but that is very noisy.

Müller et. alia define a different operator with characteristic time $\tau$ in terms of iterated EMAs $\epsilon^n_T$ of levels $n=1,2,4$ as

$$ \Delta_\tau = \gamma(\epsilon^1_{\alpha\tau} + \epsilon^2_{\alpha\tau} - 2\epsilon^4_{\alpha \beta \tau})  $$

where $ \gamma=1.22208,\beta=0.65, \alpha^{-1}=\gamma(8\beta-3)$.

This operator has a good-looking kernel with a nice small tail, but for monotonically increasing series $x$ it is not strictly positive. What differential operator works well but lacks this sign problem?  Should I just be smoothing the straight differentials?

Example

Let's say I am tracking how frequently a stock has been mentioned in the news in the last 5 years.  Some days there are no mentions, while on other days they are bountiful.  The rate of mention is of course the derivative of the number of total times the stock has been mentioned.

I could just smooth the mention counts within a window and be done, but it seems elegant to apply a differential operator to the mention count.  Unfortunately, as I noticed, I get some negative value which are clearly silly.


texts: Denoise Image with Gaussian Noise Using MATLAB / Octave

texts:
I want to remove a noise for an image using MATLAB, when the observed image is
$$f=u+v$$
where $u$ is the restored image (is the image i want recovered) and $v$ is the gaussian noise.

To restore $u$, I solve the following minimization problem:
$$\min_{u \in H^1(\Omega)} \int_\Omega \gamma|\nabla u(x)|^2dx+  \int_\Omega (f(x)-u(x))^2dx,$$
where $\gamma$ is the regularization coefficient and $\Omega=[0,n]\times[0,m]$ and $[n,m]=$size($u$).

I want to solve the PDE (Euler-Lagrange) using MATLAB:
\begin{eqnarray}
div(\gamma \nabla u) + u = f \;in \;\Omega \\
\frac{\partial u}{\partial n}=0 \; in \; \partial \Omega
\end{eqnarray}

Can anyone help me to solve this problem? Thank you!

I tried the following code :

clear all, close all,clc

uor=imread('gourd.bmp'); % the original image
u0 = imnoise(uor,'gaussian',0,0.01);
u0=double(u0);
[m n]=size(u0);
uor=double(uor);
u=u0;
c=0.028;
h=1;
for Iter=1:50,
    for i=2:m-1,
      for j=2:n-1,
          Lap=0.003*(u(i+1,j)+u(i-1,j)-4*u(i,j)+u(i,j+1)+u(i,j-1));
          u(i,j)=(u0(i,j)+(1/(2*c*h*h))*Lap);
     end
    end
    for i=2:m-1,
          u(i,1)=u(i,2);
          u(i,n)=u(i,n-1);
        end

    for j=2:n-1,
          u(1,j)=u(2,j);
          u(m,j)=u(m-1,j);
        end

        u(1,1)=u(2,2);
        u(1,n)=u(2,n-1);
        u(m,1)=u(m-1,2);
        u(m,n)=u(m-1,n-1);

en=0.0;
    for i=2:m-1,
      for j=2:n-1,
      ux=(u(i+1,j)-u(i,j))/h;
      uy=(u(i,j+1)-u(i,j))/h;
      fidelity=(u0(i,j)-u(i,j))*(u0(i,j)-u(i,j));

      en=en+c*fidelity;
      end
    end

Energy(Iter)=en;

%  Error between uor and u0
 ur=reshape(u,m*n,1);
 uori=reshape(uor,m*n,1);
 residu=norm(ur-uori)/norm(uori);

 [peaksnr, snr] = psnr(uor, u);

disp(['    iter ' num2str(Iter), ' :     Error = ' num2str(residu), ...
    ' ,    Peak-snr ' num2str(-peaksnr), ' ,    SNR ' num2str(snr)]);

  end

% show the structural similarity index for measuring image quality
[ssimval, mapssim] = ssim(u,uor);
disp([' the structural similarity index is ' num2str(ssimval)]);
figure,imshow(mapssim,[]); axis square;

figure,imagesc(u); axis image; axis off; colormap(gray);

The original image is here :
https://www.dropbox.com/s/4bccby1f4lxp4j9/gourd.rar?dl=0

Best regards


texts: remove zero padding effect crosscorrealation

texts:
I am trying to make the cross correlation on the frequency domain and I
use the fft function which takes power of 2 input data so I must pad with zero
my incoming sample signal and localy generated special code to take fft of
this signals.

I compare the first maximum and second maximum peak ratio and then I decide if the code was found or not found. When I use zero padding version this ratio and
phase change.

How can I remove zero padding unwanted effect from this correlation?

Thanks for helps.


texts: Uniform quantizer for gaussian input signal

texts:
I have a gaussian signal $x$ with zero mean and unit variance. I am quantizing it by using a uniform quantizer of step size $q$ calculated by

$$q = \frac{x_\mathrm{max} - x_\mathrm{min}}{2^{M_\mathrm{ENOB}}},$$

where $M_\mathrm{ENOB}$ is the effective number of bits. Is this the correct way of quantizing a Gaussian signal?


texts:
I'm trying to learn something more about Wiener Filters and I can't find anywhere what happens when I feed a Wiener Filter with a non WSS process. Could anyone give me an explanation or any book (specifically, which chapter/s) that explains that?

Thanks!

WSS: Wide Sense Stationary


texts: Fluctuation of autocorrelation of a signal due to signal's noise

texts:
I have a question about the fluctuation of autocorrelation of a signal due to signal's noise.

I have a signal defined in $-1\leq t \leq 1$ as the following:
$V(t)=kt+R(t)$, where $R(t)$ is the random noise, and $k$ is non-negative.

Then, when I calculate the autocorrelation of $V(t)$, $G(\tau)=\int_{-1}^{1-\tau} V(t)V(t+\tau)dt$ (using the commend CorrelationFunction of Mathematica), I could observe that the fluctuation of $G(\tau)$ reduces as $k$ increases from zero. Could you explain the origin of this observation? Is it because as $k$ becomes larger, there is more autocorrelation and $R(t)$'s contribution to $G(\tau)$ is diminished?


texts:
I have got a video containing 200 frames and i am extracting low level features at pixel level. Then i am dividing the entire frame in to small spatial patches (say of size 7x7) and creating a feature vector for each spatial patch (by summing up the feature values of the pixel inside the particular patch). If i want to use the HoG features in the above mentioned procedure, what is the best possible way to do that? i.e using PCA? or shall i use the values of the bin for that particular 7x7 patch directly (say if i have 4 bins then it will be 4 feature variables)?

PS : i am using GMM to model the features. So i want the features of each patch can be correlated with one-another using euclidean distance measure.


texts: TDE, OFDE and FDE equalizer simplest exmple in matlab

texts:
I want to implement the equalizer in matlab. I am just learning this. I want a good kick start code that can help me to learn it further. I searched a lot on Google but do not find a good reference code to start. I have a received signal in matlab that is modulated with QAM-16 modulation and passed through AWGN Noise.

I do not know much about the equalizers, I also do not know how to find impulse response of that signal. can some body please help me to just let me start some piece of code so that I can learn more in this direction.

Thanks


texts: why Sampling with higher than nyquist frequency is causing aliasing?

texts:
fm=20000 and I am sampling following signal at 65000 Hz:

x=0.5*sin(2*pi*(fm-2500)*t)+0.6*sin(2*pi*(fm-1250)*t*1.05)+0.7*sin(2*pi*fm*t*1.1)+0.8*sin(2*pi*(fm+1250)*t*1.15)+0.9*sin(2*pi*(fm+2500)*t*1.2);

and looking at frequency spectrum is not as it should look like.
Please find attachment.
It has code and result also.


texts: Measuring microphone Signal to Noise ratio (SNR) with A-weighting

texts:
Foreword: I have already read Applying A-weighting and user @endolith's  A-weighting code

I would like to measure the SNR of a microphone with the same method applied by the main manufacturers, that is displayed in their datasheets (example here with Neumann's U87 technical spec):

Here is how I plan to do it:

Play a continuous 1 kHz sine signal with a studio monitor loudspeaker (ideally in an anechoic room)
Adjust the playback volume such that the measured volume, close to the microphone, is 94 dB SPL, by using a sound level SPL meter device

Question here: should it be 94 dB SPL RMS or peak? (For a sine there's a difference of ~ 3 dB between RMS and peak)
Record the audio signal with a digital audio interface
Stop the playback of the sine and also record some silence (and thus some microphone noise: that is precisely what we want to measure!)
Normalize the recorded WAV file to 0 dB FS (so this will take the peak of the signal in consideration)

Is this step correct? It seems that a normalization is required here; if no normalization is done, the result SNR would depend on the preamp input gain!
Apply the A-weighting or CCIR-weighting (also called ITU-R 468 weighting), for example with endolith's code here
Compute the RMS of the A-weighted signal, during the silence/noise-only part.
This will be the final "Signal-to-noise ratio, A-weighted (re. 94 dB SPL) of the microphone".

Is this correct?

Here is a WAV file I recorded with the method described here. I'm curious to find how to find the A-weighted SNR of this microphone.


texts: Using superposition and time invariance to find RL circuit response

texts:
I am confusing reading book by Chaparro (Signal System with Matlab) in example 2.7. He gives a question, a RL circuit with unit-step source $ v(t) = u(t) $ will result a response,

$$ i(t) = (1- e^{-t}) u(t) $$

If the input changed $ v(t) = u(t) - u(t-2) $, we are asked to the new response of RL system.

In the given solution (see the figure), he only use superposition and time invariance principle to bring the following,

$$ i(t) − i(t − 2) = 2(1 − e ^ {−t} ) u(t) − 2(1 − e ^ {(t−2)} ) u(t − 2) $$

Based on my understanding on superposition and time invariance, it should be,

$$ i(t) − i(t − 2) = (1 − e ^ {−t} )u(t) − (1 − e ^ {(-t+2)} ) u(t − 2) $$

But, I've checked it with graphical method (convolution) and Matlab script, the given answer (plot of $i(t)$, $-i(t-2)$ and $ i(t)-i(t-2)$) is true. Here is how,

First I obtain $h(t)$ from $i(t)$ by derrivation,

$$i(t) = v(t) * h(t)$$

I got $h(t)=e^{-t}u(t)$. By convolving $h(t)$ with $i(t)-i(t-2)$ I got the same plot on the book (Using matlab conv).

Anyone can explain, how superposition and time-invariant principle can bring the answer?

Thanks in advance.

Screenshot of the book (from Google Book)

Plot of $i(t)$, $-i(t-2)$ and $ i(t)-i(t-2)$


texts: How to prove that frequency mixer is a nonlinear system?

texts: I know that frequency mixer $y(t)=x(t)cos(\omega_ct)$ is a nonlinear system since its output has frequency components that are not present in the input. To prove its nonlinearity, it shuldd be shown that the system either violates scaling or additivity. However, I could not do that. Can someone please help me?


texts: Question about vanishing moments in wavelet transforms

texts:
I am reading the book Noise reduction by wavelet thresholding by Maarten Jansen. About vanishing moments, it reads

  To create a really sparse representation, we try to make coefficients
  that live between point s of singularities as small as possible. In
  these intervals of smooth behavior , the signal can be locally well
  approximated by a polynomial.

I don't understand the sentence. If the number of vanishing moments is large, does it mean that the filtered image will be smoother? If I hope to keep edges in images, do I need to use Haar wavelets?


texts: Recognizing bodily noises (ie, not spoken words)

texts:
I was pointed here by a chain of stackoverflow posts. I'm a software developer looking to make a module that recognizes noises like snorts and coughs. Here was my original question on stackoverflow and here is where I found this SE board.

In short, I'm trying to use a Raspberry Pi and a microphone to pick up noises like snorts, coughs, sneezes, etc as a proof of concept to tell if someone is sick. I'm not quite sure where to start, as those posts I linked above have a lot of information. I've also come across this post, but I'm not sure if that helps.

I'm very new here and to this particular concept, so any and all help would be appreciated.


texts: Understanding sine wave generation in Python with linspace

texts:
I was trying to sample a 12.8 MHz sine wave (78.125 ns) signal at every 160us (micro seconds). Since 160us is multiple of base period 78.125ns(x2048) i expected to get a sample of fixed amplitude but instead what I am seeing is a another periodic sine wave. I don't understand why ?

I am doubting quantization error but shouldn't that generate uniform noise instead of a creating a periodic sine wave.

    import numpy as np
    from matplotlib import pyplot as plt

    fig2 = plt.figure()
    ax2 = fig2.add_subplot(1, 1, 1)

    capture_size1 = 2048
    timestep1 = 160e-6
    freq1 = 12.8e6
    time1 = np.linspace(0, capture_size1 * timestep1, capture_size1)
    w1 = np.sin(time1 * 2 * np.pi * freq1)
    ax2.plot(time1, w1, '.')

    plt.show()

Edit1 :
1. the 12.8 MHZ is intentionally under sampled

Adding the screenshot of the plot with capture_size1 = 2048, the sine wave has proper amplitude of [+1, -1]

Edit2: I tried to increase the precision by using Decimal and i see it is behaving as expected. I expect a straight line as the sampling point is an exact multiple of period.

from decimal import Decimal
from math import pi as mpi
from math import sin as msin
import numpy as np
from matplotlib import pyplot as plt

fig2 = plt.figure()
ax2 = fig2.add_subplot(1, 1, 1)

capture_size1 = 2048
timestep1 = 160e-6
freq1 = 12.8e6
time1 = np.linspace(0, capture_size1 * timestep1, capture_size1)
w1 = np.sin(time1 * 2 * np.pi * freq1)
ax2.plot(time1, w1, '.')

capture_size3 = Decimal(2048 * 16)
timestep3 = Decimal(160e-6)
freq3 = Decimal(12.8e6)
time3 = [Decimal(i) * timestep3 for i in range(capture_size1)]
w3 = [msin(Decimal(i) * timestep3 * Decimal(2) * Decimal(mpi) * freq3) for i in range(capture_size1)]
ax2.plot(time3, w3, '.')

plt.legend(["Actual", "Expected"])
plt.show()

Edit3:
I further did some analysis thanks to the comment by @jithin. Looks like this is an issue with linspace. I tried to generate the time interval by just multiplication as shown in below code and removed the original plot which used linspace(this is crucial) so now i am able to see the values in 1e-9 range as others suggest. So is there indeed issue with linspace ?

from decimal import Decimal
from math import pi as mpi
from math import sin as msin
import numpy as np
from matplotlib import pyplot as plt

fig2 = plt.figure()
ax2 = fig2.add_subplot(1, 1, 1)

capture_size1 = 2048
# timestep1 = 160e-6
# freq1 = 12.8e6
# time1 = np.linspace(0, capture_size1 * timestep1, capture_size1)
# w1 = np.sin(time1 * 2 * np.pi * freq1)
# ax2.plot(time1, w1, '.')

capture_size2 = 2048
timestep2 = 160e-6
freq2 = 12.8e6
time2 = [i * timestep2 for i in range(capture_size2)]
w2 = [np.sin(i * timestep2 * 2 * np.pi * freq2) for i in range(capture_size2)]
ax2.plot(time2, w2, '.')

capture_size3 = Decimal(2048)
timestep3 = Decimal(160e-6)
freq3 = Decimal(12.8e6)
time3 = [Decimal(i) * timestep3 for i in range(capture_size1)]
w3 = [msin(Decimal(i) * timestep3 * Decimal(2) * Decimal(mpi) * freq3) for i in range(capture_size1)]
ax2.plot(time3, w3, '.')

plt.legend(["multiply", "Decimal"], fontsize='xx-large')
plt.show()

The image of the above python code is below


texts: Decimating by last interpolating factor? (Multi-stage)

texts: Let's say I have a multi-stage interpolator at the transmitter, which interpolates by 2, 4 and 8 (that is, 64). When decimating (at the receiver) why should I start decimating by 8, then 4, then 2, instead of 2, 4 and 8? I mean what's the reason behind? What makes it better than doing it the other way round?


texts: Scipy.signal noise with rfft compared to fft

texts:
I'm trying to get the fourier transform of a signal with real values, however the results I get with rfft are noiser than those with fft.

I wrote the following code:

import numpy as np
from scipy.fftpack import rfft, fft
import matplotlib.pyplot as plt
# Number of sample points

N = 600
# sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
f1 = 50
f2 = 80
y = np.sin( f1* 2.0*np.pi*x) + 0.5*np.sin(f2 * 2.0*np.pi*x)
yf = fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
plt.plot(xf, 2.0/N * np.abs(yf)[0:N/2])
#plt.savefig('sin.png')
plt.show()

yf2 = rfft(y)
xf2 = np.linspace(0.0, 1.0/(2.0*T), N)
plt.plot(xf2, 2.0/N * np.abs(yf2))
plt.savefig('sin2.png')
plt.show()

and I get the following results:

I thought that for real values I would get the same result with fft and rfft, do you know why there are some differences ?


texts: sparse representation for image denoising

texts: When I read papers on image denoising, I always encounter sparse representation. For image denoising, we try to separate image signal from noise. It is assumed that signal is correlated and noise is uncorrelated. Sparse representation represents one signal as a linear combination of a small number of dictionary elements. It hope to use as few as non-zero coefficients to represent signal. In my view, it is appropriate for image compression. How can it is related to image denoising?


texts: Why PSNR decreases for Ultrasound Images?

texts:
I am working different types of filters to remove speckle noise in ultrasound image. I have taken the below paper as my base paper and trying to reproduce the results of this paper:

My Base Paper

I designed all the filters and tested with the default matlab images. I was happy with the results. The Mean Square Error (MSE) of filtered image is less than that of Noisy image and the Signal to Noise Ratio (SNR) and Peak Signal to Noise Ratio (PSNR) of the filtered image is greater than that of Noisy image.

But Sadly, When I tried it for ultrasound images the results are opposite. That is the MSE of filtered image is greater than that of Noisy Image and the SNR & PSNR of the filtered image is lesser than that of Noisy Image.

I tried for many other ultrasound images but unable to get it right.

Results for Lena.jpg

*******Quality Measurements***********

*******Noisy Image***********

Mean Square Error = 0.0080186

Signal to Noise Ratio (SNR) = 69.4875

Peak Signal-to-Noise Ratio(PSNR) = 69.0898

*******3 x 3 Median filter***********

Mean Square Error = 0.00257916

Signal to Noise Ratio (SNR) = 74.4137

Peak Signal-to-Noise Ratio(PSNR) = 74.016

*******5 x 5 Median filter***********

Mean Square Error = 0.00188434

Signal to Noise Ratio (SNR) = 75.7769

Peak Signal-to-Noise Ratio(PSNR) = 75.3792

*******7 x 7 Median filter***********

Mean Square Error = 0.00208378

Signal to Noise Ratio (SNR) = 75.34

Peak Signal-to-Noise Ratio(PSNR) = 74.9423

Results for ultrasound.jpg

*******Quality Measurements***********

*******Noisy Image***********

Mean Square Error = 0.00153502

Signal to Noise Ratio (SNR) = 64.8881

Peak Signal-to-Noise Ratio(PSNR) = 76.2697

*******3 x 3 Median filter***********

Mean Square Error = 0.00770785

Signal to Noise Ratio (SNR) = 57.8799

Peak Signal-to-Noise Ratio(PSNR) = 69.2615

*******5 x 5 Median filter***********

Mean Square Error = 0.00810142

Signal to Noise Ratio (SNR) = 57.6637

Peak Signal-to-Noise Ratio(PSNR) = 69.0452

*******7 x 7 Median filter***********

Mean Square Error = 0.00853159

Signal to Noise Ratio (SNR) = 57.439

Peak Signal-to-Noise Ratio(PSNR) = 68.8205

I am unsure why I am getting results like this. I  have added code that I have used for quality measures. Please correct me If I have done anything wrong.

function metrics = Metrics1(Orig_Image,Esti_Image)

    %---Mean-Square Error(MSE) Calculation
    Orig_Image = im2double(Orig_Image);%---Convert image to double class
    Esti_Image = im2double(Esti_Image);%---Convert image to double class
    [M N] = size(Orig_Image);%---Size of Original Image
    err = Orig_Image - Esti_Image;%---Difference between two images
    metrics.M_SE = (sum(sum(err .* err)))/(M * N);

    %---Signal-to-Noise Ratio(SNR) Calculation
    metrics.SNR = 10*log10((1/M*N)*sum(sum(Orig_Image.*Orig_Image))/(metrics.M_SE));

    %---Peak Signal-to-Noise Ratio(PSNR) Calculation
    if(metrics.M_SE > 0)
        metrics.PSNR = 10*log10(255*255/metrics.M_SE);
    else
        metrics.PSNR = 99;
    end

     %---Mean and Standard Deviation

    %---Beta Calculation
    h = fspecial('laplacian');
    I1 = imfilter(Orig_Image,h);
    I2 = imfilter(Esti_Image,h);
    I_1 = mean2(I1);
    I_2 = mean2(I2);
    metrics.Beta = sum(sum((I1 - I_1).*(I2 - I_2)))./(sqrt(sum(((I1 - I_1).^2).*((I2 - I_2).^2))));
end


texts: what does the output of ifft in ofdm look like?

texts: i have a problem in ofdm modulation .... there is a stream of bits that is converted from serial to parallel and is applied to an IFFT block. based on the definition of IFFT the output should be an complex valued number. how does it represent the real and imaginary parts in the digital domain( bits 0,1) ?


texts: Why does altering path length in radio transmission alter the phase of the received signal?

texts:
I do have a very naive yet bothering question:
Assume I send a sinusoid signal

\begin{equation}
A*cos(2\pi ft+ \phi_0)
\end{equation}

over a radio channel. Assume there is only a LOS-path for each receiver antenna and no interference or other fading effects.

Next, I receive the signal at two different locations with different distances from the receiver. Will the phase of the received signals be different for the two sites?

In order to deepen my understanding:

Let's assume I record my received signal and I set $t=0$ for each receiver individually when the signal arrives at the respective receiver.
Will the phase of the received signal at each receiver be different for different distances? Obviously, for one observation instant (time on reference clock equal for all receivers) the phases vary at the different sites. However, what I mean is:
Does the signal phase vary for the situation I explained above (individual receiver clocks)?
This question boils down to the understanding of the physical reality of an EM wave propagating. Assume I follow the wave and I 'look' at the wavefront at different times; will I see the same phase at all times? Or does the 'front' of the signal wave oscillate in time and space?
I really hope that does not sound too stupid; it's really a question that bothers me...


texts:
Am reading Rick Lyons - Understanding DSP book about Sampling. I have few doubts related to this chapter.

what is the difference between Bandpass Signal and a Passband Signal ?
Can somebody give an example for a Bandpass system which produces Bandpass signal as output !
When it comes to Baseband Sampling ( Low pass Sampling ) sampling theorem states that " Fs > twice the Highest frequency Component" and for Bandpass Sampling, sampling theorem says Fs > twice the Bandwidth of the signal. why two different sampling statements for the same signal of B hertz.

I understand in Bandpass sampling if I sample at the rate of twice the highest frequency component, then it'll be oversampling.


texts:
Given an absolutely summable signal $x[n]$, the $z$-transform $X^z(z)$ is rational with a pole at $z=0.5$.

Given the following the statements:

$x[n]$ has a finite support in the time domain.
$x[n]$ is a left sided signal.
$x[n]$ is a right sided signal.
$x[n]$ is a two sided signal.

Can anyone please tell which statements are true and explain why or why not?


texts:
All I know so far is that frequency is how much a signal repeat it self in a given amount of time (usually one second), like shown in the following figure:

I'm confused with some concepts like frequency range or I think it's the same as frequency band (correct me if I'm wrong), so what does it exactly mean?

As I understood, if the frequency range from 1 kHz to 3 kHz, does that mean that we have a signal that maybe repeat itself  from 1000 times to 3000 time every second?

I saw that the bandwidth of analog signal is the difference between highest frequency and the lowest frequency, so in the above example is 2 kHz does that mean it is the difference between  the maximum repetition and lowest repetition of the signal?


texts:
I've seen the use of a Causal wavelet mentioned a few times but have never seen any examples of them in use with Matlab.  Is it possible to to perform a Wavelet Decomposition using a causal wavelet i.e. no future data.  If so an example of how to do this would be greatly appreciated.

Thank you.


texts: Zero-pad before or after windowing for FFT

texts: What's the correct way. Should I zero-pad a signal before or after applying a windowing function?


texts: Is the modulation on my WiFi different to the modulation on cellular?

texts: Is the modulation method, or nature of modulation otherwise, between my phone and my WiFi modem, versus my phone and a cellular mast, different?


texts: Inverse Fourier Transform of piecewise function

texts:
Objective : Compute $y(t)$ from $Y(j\omega)=H(j\omega)X(j\omega)$ where :
$$
x(t)=\left(\frac{\sin(2t)}{\pi t}\right)^{2}
$$
and
$$
H(j\omega)=\begin{cases}e^{-j\omega}&\text{if $|\omega|<4$}\\0&\text{if otherwise}\end{cases}
$$
what I did is find FT of $\left(\frac{\sin(2t)}{\pi t}\right)\left(\frac{\sin(2t)}{\pi t}\right)$ as follow :
$$
\mathcal{F}\{s\cdot s\}(t)=\frac{1}{2\pi}[S*S](j\omega)
$$
We have : $$S(j\omega)=\mathcal{F}\left\{\frac{\sin(2t)}{\pi t}\right\}=\begin{cases}1&\text{if $|\omega|<2$}\\0&\text{if otherwise}\end{cases}$$
I have computed the modulation :
$$
X(j\omega)=\frac{1}{2\pi}[S(j\omega)*S(j\omega)]=\frac{1}{2\pi}\int_{-\infty}^{\infty}S(j\Omega)S(j(\omega-\Omega))\;\text{d}\Omega=\begin{cases}\frac{\omega}{2\pi}&\text{if $0<\omega<2$}\\\frac{2-\omega}{2\pi}&\text{if $2<\omega<4$}\\0&\text{if otherwise}\end{cases}
$$

The first issue I have is multiplying $H(j\omega)$ with $X(j\omega)$ since I have two piecewise functions. The second issue is how to thus compute its inverse Fourier transform.

I would hope someone can please help me. Thank you


texts:
I need some information about Stockwell transform (also known as the S-transform):

How can I implement it in MATLAB?
*Does it give you the damping ratio $\zeta$ of a signal like the complex Morlet continuous wavelet? $$x(t)= e^{-\omega_{n}\zeta}\sin(\omega_{d}t+\theta )$$
What are its advantages compared to CWT?

The S-transform defined as follows:
$$s( \tau ,f)= \int_{-\infty}^\infty h(t) \frac{\ |f|}{\ \sqrt{2 \pi } } e^{ \frac{-(t- \tau )^{2} f^{2}}{2}  }    e^{-j 2 \pi f} dt$$ (This is a good tutorial about the S-Transform)

These are the first few question that I can think about at the moment. It seems no one asked about the S-transform before, so put your questions in this post by editing it, and answer them if you know them. Someone with enough reputation put a tag for S-transform please. Thanks a lot.