texts: What does it mean to deconvolve the impulse response

texts:
I am learning DPS and I came across the problem of deconvolution and removing the impulse response from a signal? This still does little sense to me.

My understanding of the impulse response is to find out the response, usually in frequency, of a device, a filter, a room etc. If I don't mess things here, to me, the impulse response in signal processing is equivalent to numeric 1. Finding the impulse response is convolving a signal with something that is equivalent to 1 that is the impulse response. So if I need to find what the value of x is, I simply multiply x by 1:

$$y = x\times1$$

Therefore if $y = 0.8$, I know $x = 0.8$ as well.

In DSP, we have the equation

$$Y = XH$$

So to know the impulse response of $X$, I multiply $X$ by $H$ and obtain $Y$ which is equal to $X$. I've gone that far.

I don't understand what deconvolution means. I though by applying the impulse response we have already an answer about the frequency of anything we need. What does deconvolution mean here and why would we look for it?

Thanks.


texts: How to pad a windowed sinc filter in the frequency domain

texts:
I'm creating windowed sinc filters to apply them to certain signals that I'm dealing with. To design the filters, I'm using the approach described in the book "The Scientist and Engineer's Guide to DSP". Here's a brief resuming:

\[h[i] =
\left\{
\begin{matrix}
Kw(i)\frac{\sin(2\pi f_c(i - \frac{M}{2}))}{i - M/2} & 0\le i \le (M+1)\\
\\
0 & \text{otherwise}
\end{matrix}
\right.
\]

\[
w(i): \text{A window function}\\ K: \text{DC gain}\\ f_c: \text{The cut frequency expressed as a fraction of the signal sampling rate}
\]

After calculating the filter kernel, I obtain the frequency response of the filter by taking the FFT of h[i] using an algorithm that approach the filter length to the nearest power of two(to get better performance). The same procedure it's also applied to the input signal in order to get its frequency response.

To my understanding, the next step would be multiply the two frequency responses to get the filtered signal in the frequency domain, but the problem is that both signals are not the same size and I haven't found an explanation on how to deal with this problem. So people, how should I pad the filter signal to make this point-wise multiplication? Should I fill the end of the filter with zeros to match with the signal size or should I do it symmetrically putting zeros to the left and right of the filter response?

Thanks in advance!


texts: Interpreting N in DFT as the Number of Points vs. Number of Intervals

texts:
The "N" is DFT is understood to be the number of data points in a given sequence or in other words the length of the sequence. We recently have had discussions here Indexing in DFT (from an old paper) and someone's old question How do I measure the time duration of a finite-length discrete sequence?. One of the popular symbolic versions of DFT is
$$
X(j)=\sum_{k=0}^{N-1} x(k) \exp \left(-i 2 \pi\left(\frac{j}{N}\right) k\right)
$$
Suppose someone gives us the sequence consisting of N=11 points without telling us the total time or the sampling rate. If we apply DFT on it in MATLAB, the output is 11 points

Notational problems start when we wish to determine the time interval $\Delta$t and the frequency step $\frac{1}{N\Delta t}$ when the $\Delta$t is revealed.
a) If we wish to determine the sampling frequency it is (N-1) points collected in 1 seconds. The last point belongs to the next cycle. The correct sampling rate is 10 Hz not 11.
b) The frequency resolution is shown to be $\frac{1}{N\Delta t}$. In order to get the correct  frequency step, we to have to put $\frac{1}{10(0.1)}$  NOT $\frac{1}{11(0.1)}$.
Therefore, the frequency resolution, if we accept that N=11,
$$\frac{1}{(N-1)\Delta t}$$ but nobody shows this formula.
It seems that we are using N in two ways

N as the length of the sequence
N when we have to determine the frequency step where actually it is one less than N to get the correct result.

A respected colleague says that N should be interpreted as the no of intervals not the number of points. This is inconsistent with defining N as the number of points in the sequence. Unfortunately, I cannot find any reference which says that N is the number of intervals.
How can we make this consistent?
Thanks.


texts: How enhance amplitude of audio signal without distortion?

texts:
After FFT I've set to zero some frequencies and I'd like to enhance the amplitude in order to recover the leak.
I'd like to do a normalization in order to avoid distortion.
I've done in this way but the new signal is between 0 and 1 and so the result is incorrect:

maximum_data = np.max(fftdata)
minimum_data = np.min(fftdata)
maximum_abs = np.max(fftabs)
minimum_abs = np.min(fftabs)
normalization_factor = ((2**16)/2)-1
for i in range(len(fftdata)):
    fftdata[i] = ((fftdata[i] - minimum_data)/(maximum_data - minimum_data))*normalization_factor
    fftabs[i] = ((fftabs[i] - minimum_abs)/(maximum_abs - minimum_abs))*normalization_factor

How can I resolve the problem? What is the correct formula?
Thanks.


texts: Physical meaning of average values of random signals

texts:
This question might be a bit stupid, anyway, i'll risk it, since i want to get better understanding of this subject.

Let's consider random signal x(t), and let's say that we know that it is ergodic by mean, and by autocorrelation (this implies that it is wide sense stationary).

This means that time averaged mean of this signal is equal to it's statistic averaged mean, same goes for autocorrelation.

So we have  the following:

$\mu_X=E(x)=\frac{1}{T} \int_T x(t)dt=\langle x(t)\rangle  \\ R_{xx}(\tau)=E[x(t)x(t+\tau)]= \lim_{T->\infty} \frac{1}{T} \int_T x(t)x(t+\tau)dt=\langle x(t)x(t+\tau)\rangle$

Then, as i found in one of the books i have, it says following:

$\langle x(t) \rangle$- corresponds to DC level of given signal

$\langle x(t) \rangle ^2$ - corresponds to normalized power of DC component

$\langle x^2(t) \rangle$ - corresponds to total average normalized power

$\overline{\sigma_x^2}= \langle x^2(t) \rangle - \langle x(t) \rangle^2$-corresponds to average normalized power in AC component of the signal

$\overline{\sigma_x}$-corresponds to rms of AC component

My question is, is there any way to prove this, or to gain intuitive understanding why this is true. Any help appreciated!


texts: Deriving the transfer functions of a heating system

texts:
I'm tying to develop a model for a heated system that consists of a small steel block (~25 sq. in.) with a heating element embedded in it. The block acts as a sort of hot-plate that is used to deform other materials. There is a temperature sensor embedded in the block as well for feedback. I'm only trying to model the temperature of the block, not the temperature of the deformed materials.

I've used Newton's Law of Cooling to determine the time constant for the steel block cooling with open air convection:

$$
{dT(t) \over dt} = -k \left(T(t) - T_{\rm ambient}\right)
$$

where $k$ is the time constant (and is ${\approx 0.00228} sec^{-1}$). Taking the Laplace Transform of this gives me:

$$
H(s) = {T(0) \cdot s - k \cdot T_{\rm ambient} \over s - k}
$$

I've also been able to determine the heating capacity of the heating element coupled with the steel block:

$$
{dT(t) \over dt} = -k \left(T(t) - T_{\rm ambient}\right) + Q
$$

where $Q$ is the heating (and is ${\approx 0.412}$ $°C/sec^{-1}$). Taking the Laplace Transform of this gives me:

$$
H(s) = {T(0) \cdot s - k \cdot T_{\rm ambient} + Q \over s - k}
$$

I've simulated the responses of both of these transfer functions under step inputs, and they match closely with experimental results. But when I attempt to add closed loop control, I get simulation results that do not make sense to me. The transfer function I've been calculated for simple negative feedback is:

$$
H(s) = {{T(0) \cdot s - k \cdot T_{\rm ambient} + Q} \over {(T(0)+1) \cdot s - k \cdot T_{\rm ambient} + Q - k}}
$$

Simulating this transfer function with a step input produces a system that goes to the target temperature almost immediately. I'm expecting to see an exponential response starting from the $T(0)$ point and increasing to the target temperature, similar to what I see when simulating the open loop responses.

I've gone back to my old textbooks and searched the internet for additional resources. I've solved a number of mass-spring and RC circuit problems, but I can't seem to make the connection between those problems and this one. So I'm not sure where I'm going wrong.


texts: Finding Z transform of a signal: Intermediate steps

texts:
Find the Z transform of

  $y(n)=x(n+2)u(n)$

I have solved the problem. I have doubt whether it is correct or not. It would be very helpful if someone could check whether the steps that I have done are correct or not. Thank you.


texts:
Which method provides better noise immunity?

I know you can transmit more data via psk, but what about which method is less susceptible to noise?


texts:
I have two variables of data (BxHPF and ByHPF), which are equall in number of data. I would like to make a loop, which takes first value of BxHPF and ByHPF compute them, then take next one from BxHPF and ByHPF and so on. For now I have:

So I though about this way:

Computing = (float(BxHPF[0])/float(ByHPF[0]))
Dat1 = math.degrees(Computing)*(-1)

In fact those equations give me correct result. But as I wrote I need to make a loop, which will count every pair from BxHPF and ByHPF with using those two variables Computing and Dat1.

To be exact BxHPF and ByHPF contain 266150 records each.


texts: Why does FEC seem to make BER worse at low Eb/N0?

texts:
I've been plotting BER curves and am wondering why FEC seems to make a BER curve worse at low Eb/N0 than the un-coded theoretical BER curve?

My example plot below is for Differential BPSK using a Rate 1/2 k=7 convolutional code on the AWGN channel.


texts: Phase noise estimation and compensation schemes

texts:
Brief background

There is a known problem that high order modulation schemes suffers from phase noise of reference generator. It's very important issue for e.g. line-of-sight (LOS) modems exploiting QAM-256 or even higher. Some compensation schemes are possible to solve this task. The main challenge in these schemes is to estimate phase noise in a time span of data reception. Knowing actual samples of phase noise we can eliminate it from data symbols significantly improving BER for high order modulation. Such estimation is possible due to phase noise is in general low band random process which can be adaptively predicted/estimated. One can trade off communication speed inserting known pilots in data flow to estimate phase noise at some moments of time. Phase noise from pilots' time span could be used to evaluate phase noise for the interval of data samples. So the question is about pilot mapping scheme and estimation techniques.

Point

Gathering some information about this topic, I can pick out two major approaches:

1) Sending a group of pilots (QPSK) $\to$ eliminating pilot modulation to get phase noise (PN) samples $\to$ learning of adaptive prediction filter (e.g. by LMS) $\to$ finding predictive PN samples $\to$ eliminating PN from data samples $\to$ wait for the next pilot group and repeat. It's extrapolation method.

2) Inserting pilot sample (QPSK) periodically in the data stream, e.g. 1 pilot/32 data $\to$ evaluating PN at pilots (like discussed above) $\to$ finding PN between pilots by interpolation, number of pilots from left and right are used $\to$ eliminating PN from data samples. It's interpolation method.

Has someone any experience with similar technologies? Or maybe someone face to phase noise compensation problem in practice?

I wonder to know which of these approaches (I mean pilot-based, interpolation or extrapolation) is more preferable/robust in real life modem design?

Remark

Work environment is to be considered as a channel with AWGN. SNR is pretty high, e.g. 40-50 dB. Phase noise is in the order of -70...-90 dB/Hz. Phase noise bandwidth is significantly less then signal's one.

Here is an interesting article in IEEE:

  Simon V., et el. - Phase noise estimation via adapted
  interpolation

It's about the second approach. Adaptive interpolation via LMS is discussed. Does someone know any interesting references, I'm curious to learn about.

In the picture above the first case is illustrated.

Thanks in advance


texts: How many things can they be done with a spectrogram?

texts:
I found some things that can be done with a spectrogram.

filter frequencies by setting the bins to zero
observe what frequencies make up the signal.
observe the energy or
amplitude of each frequency, the whitest pixel has more volume and the black pixels have little amplitude.
observe harmonics

It does not occur to me that other things can be done.


texts: How many non-zeros samples in convolution

texts:
If I have $b[n]$ of length 30 (30 non-zero samples) and $c[n]$ of length 40 (40 non-zero samples).

How many non-zero samples will $a[n]=b[n] * c[n]$ have? (Note '$*$' is a convolution).

I think there will be 30 non-zero samples. Am I correct?


texts: How channel state information is calculated from Sounding Packet

texts:
What I understand is for Explicit Transmit Beamforming scenario, the Transmitter sends an NDP packet(aka HT Sounding) which consist of OFDM training symbols in packet preamble (HT-LTF) which are randomly generated +/-1 for each subcarrier(say 56).

I do not have much knowledge on this, but if anyone can help me to understand how these training symbols are calculated at the receiver end to evaluate CSI will be helpful.


texts: signal analysis and synthesis with different window length and FFT length

texts:
I'm coding an analysis/synthesis system on MATLAB. I was able to recover the original signal when the length of the window (winLen) is set equal to the length of the FFT (FFTLen). When they are not equal, the results becomes erroneous. What is the cause of this error?

The analysis side consists of a time-domain input which is windowed with a Bartlett window followed by a FFT. The window shift corresponds to half of the window length. The FFT-ed sequence is decomposed into log magnitude and phase components. (Since the sequence is real, only the first half of the magnitude is stored.)

The synthesis side consists of concatenating the magnitude with its flipped version multiplied by the phase, which gives the initial FFT results before decomposition. Then IFFT is performed where the length is equal to the length of the time-domain window. The result is obtained from OLA process (since the Bartlett window is used in the analysis stage, a simple OLA will yield the original sequence).

% load input
cc = cc(1,:);
[~,T] = size(cc);

% window parameters
winLen = 32; % perfect reconstruction when winLen=FFTLen
winShift = winLen/2;
FFTLen = 64;
N = ceil((T-winLen)/winShift) + 1; % number of frames

% ANALYSIS PART
ms = zeros(1,(FFTLen/2)+1,N);
mp = zeros(1,FFTLen,N);
y = [cc(1,:) zeros(1,(winLen + winShift * N) - T)]; % zeropad sequence
for n = 1:N
    ibeg = 1 + (n-1) * winShift;
    iend = winLen + (n-1) * winShift;

    buffer = y(ibeg:iend) .* bartlett(winLen)';
    buffer = fft(buffer,FFTLen);

    ms(1,:,n) = log(abs(buffer(1:(end/2)+1)));
    mp(1,:,n) = unwrap(angle(buffer));
end

% SYNTHESIS PART
ccinv = zeros(1,winLen + (N-1) * winShift);
for n = 1:N
    ibeg = 1 + (n-1) * winShift;
    iend = winLen + (n-1) * winShift;

    msbuf = exp([ms(1,:,n) fliplr(ms(1,2:end-1,n))]) .* exp(1i * mp(1,:,n));
    buffer = ifft(msbuf,winLen);

    ccinv(1,ibeg:iend) = ccinv(1,ibeg:iend) + buffer;
end
ccinv = real(ccinv);

% PLOT RESULTS
figure(1)
index = 50:100;

hold on
plot(cc(index),'k:')
plot(ccinv(index),'r:')


texts: Discussion - Integrating Accelerometer data to Position data from Frequency domain FFT

texts:
I came across one blog who suggested to convert from accelerometer data to position data from frequency domain. Here are the steps that the original individual suggested:

  Remove the mean from the accel. data
  Take the Fourier transform (FFT) of the accel. data.
  Convert the transformed accel. data to displacement data by dividing each element by -omega^2, where omega is the frequency band.
  Now take the inverse FFT to get back to the time-domain and scale your result.

I am wondering if anyone have ever used this method before, and may I know your algorithm? I am not particularly sure on what does it means by Step 3 - frequency band?

p.s. I hope this question is relevant here.


texts: THD: estimation using periodogram or fft?

texts:
I want to estimate the THD of a voltage signal (output from a frequency inverter). At the begining, I need to estimate the power of the signal at frequency domain. To doing so, I should use a periodogram or a FFT.

Matlab uses periodogram in his THD function. I know that are differences between FFT and periodogram, as shown in this question. But, wich method should I use? and why?

Thanks.


texts: Solving Wiener Hopf integral equation for causal filter of predictor

texts:
Given a stochastic signal $x(t)$ with autocorrelation function $R_{xx}(\tau)=\mathrm{exp}(-
\alpha|\tau|)$, $\alpha>0$.

I want to predict $x(t+\lambda)$,$\lambda>0$  by $x(t-\tau)$, $\tau\ge0$ through a linear etimation,

i.e. I want to find estimator $\hat{x}(t+\lambda)=\int_0^{\infty}h(x)x(t-x)\mathrm{d}x$.

By the principle of linear minimal mean squared error, I can get Wiener Hopf equation of the integral form:
$$R_{xx}(\tau+\lambda)=\int^{\infty}_{0}h(x)R_{xx}(\tau-x)\mathrm{d}x $$
$$\implies \mathrm{exp}(-\alpha|\tau+\lambda|)=\int^{\infty}_{0}h(x)\mathrm{exp}(-\alpha|\tau-x|)\mathrm{d}x$$
I tried to use Fourier transform to solve this integral equation for $h(x)$

$$\frac{e^{j\omega\lambda}}{\\\alpha^2+\omega^2}=H(\omega)\times \frac{1}{\\\alpha+j\omega} $$

But I am not sure the right-hand side is correct or not. Is it legitimate to view right hand side as convoution of $h(\tau)$ and $\mathrm{exp}(-\alpha\tau)\mathrm{u}(\tau)$, where $\mathrm{u(\tau)}$ is a unit step function?
The textbook by Papoulis and Garcia use the factorization teqnique (it's on p594 of Papoulis 4th ed and p616 of Garcia 3rd ed) but I cannot see how to factor the transfer function to causal and non-causal part. I guess the answer of $h(x)$ is $c\times\delta(x)$ for some constant. Moreover, I want to find the mean sqaure error of this predictor.


texts: At which point in the signal processing chain to use Hann Window?

texts:
I have a question about signal processing. Experienced with embedded design and code, but DSP is in new territory.

Thus far I am using the CMSIS DSP library functions, which are working well. I am sampling a low frequency signal. Using a timer interrupt, the sampling rate is 40x higher than the highest signal I expect to see. The "buffer" size I process is 2048. I first apply a low pass FIR (arm_fir_f32), and then to the arm_cfft_f32 function so I can then calculate the frequency information. It works well enough to a point where the frequency is +/- 1Hz. However, I am hoping to get this to +/-200mHz.

I understand that I need to reduce the power (noise) at the endpoints, and the Hann window should be the best (better?) choice. Thus my questions:

At which point do I apply the Hann algorithm? (as follows)

$$w[n] = 0.5 - 0.5 \cos \left( \frac{2 \pi n}{N - 1}\right)$$

Should $N$ be the size of my sample block?

EDIT:

Code for windowing function:

(INPUT_SAMPLES=2048)
for(uint16_t i = 0; i < INPUT_SAMPLES; i++)
{
    ADC_Values[i] = 0.5 - (0.5 * cos ( (2.0 * PI * ADC_Values[i]) / (INPUT_SAMPLES - 1)));
}

Thank you for any help on this,
Roger

EDIT 2:

jojek,

creating the window and changing the calculation as suggested has fixed it. Although the window has not improved resolution as much as I'd like, it's better, and close enough to be quite usable. Fantastic...thanks for all the replies everyone.


texts: Is there a rule of thumb for zero padding in image processing?

texts:
I see there are a lot of answers on why zero padding is necessary and how it avoids wrapping around the sides of images. However is there a rule of thumb on how much padding will be good for the image processing? I am mainly thinking from the point of view of:

a.) speed of processing the image

b.) noise

I have till now referred to the following answers:

Why images need to be padded before filtering in frequency domain
Advantages/disadvantage of zero padding
FFT Zero Padding - Amplitude Change

Response will be really appreciated.

Thank you.


texts: How should I select window size and overlap?

texts:
I'm new to DSP and I have a wav file with a number of different frequencies playing at different volumes. I'd like to plot the power of the 350Hz frequency signal at each point in time. The source is rotating so should see the power rising and falling periodically. To start with I broke the signal up into windows and then implemented Geortzel's algorithm (from Wikipedia). This seems to work: https://www.dropbox.com/s/aq1aigdx8dq306q/Screenshot%202015-03-18%2017.38.20.png?dl=0 but for different window length and window overlaps I get different plots.

I'd love some guidance on what the correct way to choose the window length and the amount by which each window overlap is? window_length and window_skip in the code.

Also should the windows be processed (i.e. convolved with a window function like Hanning) I've read about this but I'm not clear on what this would do to the output, nor how to choose the length of the window function?

Here's the audio file: https://www.dropbox.com/s/knvwqwg7s3mcz56/example.wav?dl=0

Here's my code, comments welcome:

import pylab
import scipy.io.wavfile
import numpy as np
import math

def goetrzel(data, target_frequency):
        s_prev = 0
        s_prev2 = 0
        normalized_frequency = 2.0 * np.pi * target_frequency / len(data)
        coeff = 2.0 * np.cos(normalized_frequency)
        for sample in data:
                s = sample + coeff * s_prev - s_prev2
                s_prev2 = s_prev
                s_prev = s
                power = s_prev2 * s_prev2 + s_prev * s_prev - coeff * s_prev * s_prev2
        return power

def sliding_window(data, length, skip=None):
        if skip is None:
                skip = length
        n = len(data)
        for i in range(0, n - int(length), int(skip)):
                yield data[i:i+length]

detect_frequency = 350.0
#detect_frequency = 485.0
#detect_frequency = 633.0
#detect_frequency = 796.0

rate, data = scipy.io.wavfile.read('example.wav')

window_length = 2 * detect_frequency # How should I pick this?
window_skip = window_length # How should I pick this?

power = [goetrzel(window, detect_frequency) for window in sliding_window(data, window_length, window_skip)]
pylab.figure()
pylab.plot(power)
pylab.show()


texts:
I'm reading this paper on a 3-phase rectifier control.  3-phase rectifiers are usually connected to 50-Hz or 60-Hz grid. One popular control strategy is to transform the 3-phase currents and voltages to the DQ axes. Basically, the 3 voltages are changed in two "DC" voltages named $d$ and $q$. Same thing for currents, the 3 currents are changed to $d$ and $q$.

As the $d$ and $q$ values are now "dc" values. One can use PI controllers in order to reach a steady-state error of 0. So far so good. However, keep in mind that in this new axis, a frequency of 0 Hz corresponds to a frequency of $\pm 60$ Hz in the real world. (or $\pm 50$ Hz)

Here's my problem :

They manipulate the transfer functions as if they were on the same basis

For instance the Process transfer function is :

$$G_p(s) = \frac{1}{R+sL}$$

The process transfer function is in the "real-world" where DC is 0 Hz.

Then they apply an internal-model control to design the controller

Controller : $$G_c(s) = Kp+\frac{Ki}{s}$$

They use $$Kp = \frac{L}{\lambda}$$ and  $$Ki = \frac{R}{\lambda} $$

However, the controller works on the $dq$ axes where 0 Hz is $\pm 60$ Hz.

Am I wrong in thinking that you cannot simply design the controller this way?

Edit : Just to clarify, I don't doubt that their gains are good and yield a stable controller. I just doubt the rationale behind the gain calculations.


texts: How to match plot with formula in DNL INL

texts:
by definition DNK(K)=(V_out(K+1)-V_out(K)-Vlsb)/Vlsb
So for K=2
DNL(K)=(V_out(K+1)-V_out(K)-Vlsb)/Vlsb=((3/8)-(3/8)-(1/8))/(1/8)=-1 [LSB]
INK(K)=(V_out(K)-KVlsb)/Vlsb=(3/8)-2(1/8))/(1/8)=1
In the place Where DNL is +1 i get -1 for some reason.I dont know how to interpet this large jump form DNL INL.
Howdo i methamichly derive the INL DNL of these K=2 case?
Thanks.


texts:
What is the definition of Generalized Cross Correlation?

I need it for estimating the delay (time offset) between two audio signals. The time offset between the two signals might be as large as 500ms.


texts: How to choose suitable segmentation process for image processing

texts:
Hello I'm new on image processing.

I know there are a lot of threshold and region growing algorithms to use.

But I would like to know what is the trick for choosing suitable method for an image ?

Currently I'm trying to extract characters from this image

and most of the methods I've tried didn't really satisfy me. So I thought maybe I am doing it wrong. Here are some samples.

Normal Threshold

Sobel

Adaptive Threshold


texts: The Standard Deviation of The Derivative of a Signal

texts:
Given a signal with zero mean and a standard deviation of 0.1 sampled at 5000 Hz.
What would be the Standard Deviation of its 1st, 2nd and 'n' derivative?

For instance, let's say we measure the distance to a constant range target.
The STD of the sensor measurements of the range is 0.1 [Meter] (Assume the error is white) and its sampling rate is 5000 Hz.

What would be the standard deviation of the acceleration process?

Thanks.


texts: How doesn't JPEG huffman table contain 256 codes?

texts:
Here are the DHT entries I got from JPEG snoop:

  Destination ID = 1
  Class = 1 (AC Table)
    Codes of length 01 bits (000 total):
    Codes of length 02 bits (002 total): 00 01
    Codes of length 03 bits (001 total): 02
    Codes of length 04 bits (002 total): 03 11
    Codes of length 05 bits (004 total): 04 05 21 31
    Codes of length 06 bits (004 total): 06 12 41 51
    Codes of length 07 bits (003 total): 07 61 71
    Codes of length 08 bits (004 total): 13 22 32 81
    Codes of length 09 bits (007 total): 08 14 42 91 A1 B1 C1
    Codes of length 10 bits (005 total): 09 23 33 52 F0
    Codes of length 11 bits (004 total): 15 62 72 D1
    Codes of length 12 bits (004 total): 0A 16 24 34
    Codes of length 13 bits (000 total):
    Codes of length 14 bits (001 total): E1
    Codes of length 15 bits (002 total): 25 F1
    Codes of length 16 bits (119 total): 17 18 19 1A 26 27 28 29 2A 35 36 37 38 39 3A 43
                                         44 45 46 47 48 49 4A 53 54 55 56 57 58 59 5A 63
                                         64 65 66 67 68 69 6A 73 74 75 76 77 78 79 7A 82
                                         83 84 85 86 87 88 89 8A 92 93 94 95 96 97 98 99
                                         9A A2 A3 A4 A5 A6 A7 A8 A9 AA B2 B3 B4 B5 B6 B7
                                         B8 B9 BA C2 C3 C4 C5 C6 C7 C8 C9 CA D2 D3 D4 D5
                                         D6 D7 D8 D9 DA E2 E3 E4 E5 E6 E7 E8 E9 EA F2 F3
                                         F4 F5 F6 F7 F8 F9 FA
Total number of codes: 162

And

  Destination ID = 1
  Class = 0 (DC / Lossless Table)
    Codes of length 01 bits (000 total):
    Codes of length 02 bits (003 total): 00 01 02
    Codes of length 03 bits (001 total): 03
    Codes of length 04 bits (001 total): 04
    Codes of length 05 bits (001 total): 05
    Codes of length 06 bits (001 total): 06
    Codes of length 07 bits (001 total): 07
    Codes of length 08 bits (001 total): 08
    Codes of length 09 bits (001 total): 09
    Codes of length 10 bits (001 total): 0A
    Codes of length 11 bits (001 total): 0B
    Codes of length 12 bits (000 total):
    Codes of length 13 bits (000 total):
    Codes of length 14 bits (000 total):
    Codes of length 15 bits (000 total):
    Codes of length 16 bits (000 total):
    Total number of codes: 012

For AC table we are encoding the RRRRSSSS byte that contains zero run length and magnitude of nonzero AC coefficient that follows. This byte has 256 possible values. Then for DC we are encoding the SSSS nibble that can has 16 possible values.

Why then is it that in all the JPEG images I have seen, downloaded or created in MS Paint, the AC table has only 162 entries and DC has only 12? This means that for certain values of RRRRSSSS and SSSS, there is no way to compress them.


texts: How to simplify the Fourier Series Using an Approximation?

texts:
I have a signal, $f(t)$.  I know a function that can be used to generate this signal, such that I can determine its Fourier series.  I want to express this Fourier series in simpler terms so that the waveform can be rendered by a GPU in parallel, with the added ability to change the lowpass cutoff frequency of the signal dynamically.

Take for example a saw wave with a fundamental frequency of 1 Hz, represented by the function:

$$f(t) = {-2\tan^{-1}(\cot(\pi t)) \over \pi}$$

I'm using this form of the saw wave because it's cyclic and easy for a GPU to calculate.  Its Fourier series is rather difficult to solve, but I have the series on hand:

$$g(t, n) = {-2 \over \pi} \sum_{k=1}^n {\sin(k \pi t) \over k}$$

Now, I could technically use additive synthesis to generate the waveform using this series, but my project has real-time constraints and this simply won't do.  A single saw wave might have hundreds of partials up to the Nyquist frequency.  This takes a lot of time to compute, and on a GPU it's a branching nightmare.  Since I'm rendering all of the samples in parallel, there is no history buffer for existing filtering algorithms to operate on, and I can't do that either.

What I'm trying to do is approximate this series as a function of $n \in Z$, without the summation.  That way, it doesn't matter what cut-off I use, I get consistent performance across the board that's most likely a huge improvement over additive synthesis.  Right now I'm trying to fit the curve manually, but this is a laborious process.  Is there a way to find $h(t, n) \approx g(t, n)$, where $h(t, n)$ is the approximation of the series?


texts: Changing voice characteristics so it sounds like someone else voice

texts:
I have a recording of speech, and I want to manipulate the voice, so a listener
will perceive the manipulated voice and the original one as belonging to two
different people. The manipulated voice should sound very natural, not computer-generated.

Just to clarify, the manipulated voice does not have to be modeled after the voice of someone; I'm not asking about mimicking a voice, as discussed here.

I have already used the SoundTouch library for pitch shifting, and it produces very nice results, but I want to change not only the pitch.

I saw this impressive video of a software called Voice Sculptor (couldn't find the company site), which demonstrate such a manipulation.

So I'm asking on a conceptual/algorithmic level how this can be done
And also if you know how to achieve this with specific tools.

Last note: I'm new to DSP, so I will appreciate any explanation of why and how.


texts: Is MATLAB the best coding tool for Digital Image Processing?

texts:
I am a beginner. I do coding in Matlab only for image processing purpose. But my friends say it is very slow and in OpenCv I could get multiple times faster output results compared to Matlab.In basic reference book like DIP by Gonzalez and Woods ,they have used Matlab for programming.

So is Matlab the best coding tool for Digital Image Processing?
What are the advantages and disadvantages of coding in  Matlab? Also,What are other good options by which we can do coding for image processing ?


texts: Playing sounds with only one frequency at a time

texts: We have a "speaker", that can make a single frequency at a time (0 - 130k Hz). So we can play a tune by making it play a tone for 1 sec and then the next tone for 0.5 sec etc. We can change the length and frequency of the tones pretty quickly. But we can play any sound as they are typically a superposition of various frequencies (Fourier)... Or can we? Is there some algorithm to convert a sound signal into a sequence of sine waves.


texts: Smoothing 3D Data for the Second Derivative

texts:
I am not a signal processing expert and my feeble attempts at solving this problem have come up short.

I have a C++ application which is being fed regularly-spaced (in time) 3D position samples.  The samples are generally accurate and, if plotted, appear to be quite smooth.  However, when the second derivative is computed (to look at accelerations), there are occasional/rare spikes in the data which are unreasonable.  (I know the general range of possible accelerations.  And by "rare", I mean that over about two ours of data, I generally see less than a dozen spikes which are out of the realm of possibility.)  I want to do as little to the data as possible, but smooth it enough to remove the spikes.

My first attempt was to use a Kalman filter in our math library.  This removed the spikes, but changed the data far more than what I would really want and seems heavy-handed for what I need.

There seem to be a lot of options out there for 2D data, but I don't have enough signal processing background to understand the trade-offs or necessarily what might translate well to 3D data.  Any education or suggestions would be appreciated!


texts: Difference between MATLAB's $\tt butter$ and my calculated coefficient

texts:
So I figured recently that MATLAB's butter function was giving slightly different coefficients than those I calculate for a simple first order low pass filter, and I am not so sure why.

Here's how I calculate, starting from the following definition of a simple low passe filter and Bilinear Transform.

$$G(s)=\frac{1}{1+\frac{s}{\Omega_c}}$$
After bilinear transform, I get

$$ H(z) = G\bigg(\frac{2}{T}\frac{z-1}{z+1}\bigg) = \frac{1}{1+\frac{\frac{2}{T}\frac{z-1}{z+1}}{\Omega_c}} = \frac{\Omega_c}{\Omega_c+\frac{2}{T}\frac{z-1}{z+1}}$$
Then
$$ H(z) = \frac{T\Omega_c(z+1)}{T\Omega_c(z+1)+2(z-1)} = \frac{T\Omega_c(z+1)}{T\Omega_cz+T\Omega_c+2z-2} = \frac{T\Omega_c(z+1)}{(T\Omega_c+2)z+(T\Omega_c-2)}$$
From there I define a substitution variable $k=\frac{T\Omega_c}{2}$
$$H(z) = \frac{2k(z+1)}{(2k+2)z+(2k-2)} = \frac{k(1+z^{-1})}{(k+1)+(k-1)z^{-1}}$$

Now, I compare my results with MATLAB's butter like this:

function fof( fc,fs)
    wc = 2*pi*fc;
    k = wc/(2*fs);
    b = [k k];
    a = [(k+1) (k-1)];
    b = b./a(1);    % Normalize coefficients for comparison only.
    a = a./a(1);

    [b2,a2] = butter(1,2*fc/fs);
    freqz(b,a,10*fs,fs);
    hold on
    freqz(b2,a2,10*fs,fs);
end

The bigger the cutoff frequency, the bigger the difference. It looks like MATLAB is indeed right as the freqz always shows a $-3\textrm{ dB}$ gain for the butter filter, and not mine.

Where am I wrong ?


texts:
Given a captured noisy BPSK signal, how do I formulate the maximum-likelihood problem so as to estimate the matched filter? The model of interest would be baseband BPSK with AWGN:

$$c(t) = \sum_{\ell=-\infty}^\infty a_\ell s(t - \ell T) + N(\mu,\sigma)$$

Where $a_\ell \in \{-1,1\}$ and $s(t)$ is the pulse function that I would like to estimate from samples of $c(t)$.

There are a lot of maximum-likelihood estimators for the different signal parameters (delay, phase, and the symbols themselves), but I've not found ML treatment of the pulse shape itself.


texts: can we use machine learning concepts to extract blood vessel from retinal image?

texts: I am working on blood vessel extraction, i want to know whether we can extract blood vessel using machine learning algorithm? which algorithm i can use to yield good results?


texts: How to compute inner product of Wavelet transform convoluted with signal

texts:
I have two datasets $X_1$ and $X_2$ in a sparse wavelet basis, and I have two filters $f_1$ and $f_2$. I’d like to compute the inner product of the convolutions $$\langle X_1 \star f_1, X_2 \star f_2\rangle$$ in the sparse wavelet basis. How do I this?
Thanks!
James


texts: Correction of a signal through a transmitter

texts:
I am inquiring as to a practical way to solve a problem I have.

Basically, I need to transmit a signal, $x[n]$, through a seismic transmitter. (It will go through a D/A, etc). The transmitter that we have has a certain frequency characteristic, in that, it passes certain frequencies out better than others. (Nothing shocking, but for the sake of simplicity we can assume it is linear across frequency).

My signal $x[n]$ is somewhat wide-band. Hence, at the output of the transmitter, it might become distorted.

I am wondering, knowing my ideal signal $x[n]$ and knowing my transmitter frequency characteristics, is there a way I can 'massage' $x[n]$ before putting it through my transmitter, such that at the output I get the original $x[n]$ that I wanted?

Like I mentioned I am seeking guidance on practical and actionable techniques for something like this.

If the context helps this is for a seismic transponder app.

Thanks!


texts: Accelerometers (again...) using appropriate transformation for drifting signal returning back to x-axis (0 m/s)

texts:
Quick question on processing velocity results from accelerometer data. I'm aware of the unreliability of the data, but I'm willing to tolerate a fair amount of error. Also willing to perform lots of corrections...

I currently have a velocity graph from integrating acceleration which looks like as below:

And the true motion of the physical event looks more like:

So it is very interesting that the accelerometer does contain the correct signal, just that, notoriously, the accelerometer has no reliable coordinate frame and appears to lose orientation of the movements.

So my question is, to the knowledgeable signal processors on here, is there a simple correction function which would perform this transformation to get this curve to flatten down to 0 m/s?

I have fit linear lines & subtracted those from the data to fix the start of the event to the x-axis... I know I could pick the point where the data begins to drift and fit another linear line and subtract that error again. But, surely there is something more elegant?

(Other prior graphs which might be relevant:

You can see the drift begin to happen in the acceleration data basically right in the middle of the movement, the offset starts. Otherwise, I would assume, the acceleration data would be a complete mirror if bi-sected...)


texts:
Faster than Nyquist signaling is used to improve the spectral efficiency by reducing the time spacing (relaxing the orthogonality constraint) to pack more data in the same channel while tolerating a certain intersymbol interference (ISI).

What is exactly faster than Nyquist? Is it by increasing the sampling rate higher than the Nyquist rate? But this doesn't create ISI. What is the relationship between the Nyquist rate and the spacing in the time domain?


texts: How to get the background noise of microphone?

texts: I am working on the automatic volume control for television which should take care of background noise. For this, I need to first find out the background noise. I have read in research papers that we need to subtract original audio signal which can be obtained from audio amplifier present in tv from the recorded signal by the microphone. Here microphone will be recording both original audio signal plus noise. But they have also mentioned about phase and amplitude correlation. I do not understand how to do that.


texts:
Currently our software is detecting tones using FFT and checking if there is only a single peak in frequency spectrum. Our software needs to know what kind of input it is before taking action. So we differentiate the case when tone is fed as input or if the input is music.

However we are trying to reduce the MIPS and wondering if there is any other better ways to detect tones.

I found correlation with previous data can be used but again don't know the right method to do it. Is there any better method?


texts: Is an extra symbol at the origin of a QAM constellation ok?

texts:
This question is specifically directed at Maximilian Matthé (as you work with 5G systems), and to anyone who can shed some light.

I am trying to implement a GFDM system with a total of N = K*M symbols, where K = total number of carriers, and M = total number of subsymbols. I believe that in order to improve the Out Of Band radiation only a subset of carriers, Kon < K, are allocated to carry information while the other carriers are left as blanks.

Therefore, when mapping the data symbols (that are QAM modulated) of the allocated carriers into a 2D time frequency matrix of the size KxM, there are zeros introduced into the modulated data. This results in an extra symbol at the origin of the QAM constellation.

This seems incorrect as there is another symbol that is part of the constellation. For example, 16-QAM would now have an additional symbol at the origin as shown in the image below.

However, I feel like this might be ok because the extra symbol is at the origin, and therefore has no amplitude, phase, and surely would not require any transmission energy.

Although, on the receiver side the constellation looks like this:

So now the added symbols at zero seems to have some amplitude and phase. Would they be taking up some energy ?

So my question is, is this ok? It is justifiable at all?
Or should i rather skip out adding these zeros ?


texts: Help with filter stability, ROC and output of a system

texts:
I am trying to solve these two questions, but I am stuck with the filter stability and ROC. Hope somebody could help me with it.

According to what I understand, the third option is correct. Fourth option is also incorrect as it only contains one pole and on zero. I think the first and second options are incorrect; the system is unstable and ROC doesn't include unit circle because the pole lies at z = -2 and the system is casual. Not sure though. I'd appreciate some help in this regard or if any of the other options are incorrect.

For this question, the pole is outside the unit circle, so I thought the system must be unstable, but then I realize the output depends on future input (x[n+1]), which makes it non-casual, also number of zeros>poles. Since for a non-casual, the ROC moves inward, this system includes unit circle and is stable? But doubtful about the second part of the option regarding being dependent on input.
I wrote the following code to find y[n],
n = 100
w = np.linspace(-np.pi,np.pi,n)
z = np.exp(1j * w)
H = z* (z - 0.5)/(z + 2)
t = np.arange(-n/2, n/2)
x  = sp.unit_impulse(n) + 0.5 * sp.unit_impulse(n,1)
plt.plot(t,ifftshift(ifft( H * fft(x))))

Got the following result,

I don't think, 1st, 2nd and 5th options are correct, looking at figure. As I increased the value of n, the peak became narrower, so for the 4th option, I don't know if it is showing diverging oscillation. The only peak in this question is at zero.
Any help is appreciated.
Edit: Question 1 seems to work with the above logic. Stuck with no. 5 now.


texts:
What is the ordering of size  following image types : RGB,grayscale, indexed image, binary image?

RGB highest size?and binary lowest?


texts: Determine causality and stability from given filter structure

texts:
I have the diagram above. I found the transfer function below from it;

The question asks me to find out if the system is causal and stable, but didn't it have to indicate whether it was left-sided or right-sided?


texts: How to measure noise, consonance, dissonance of a sound?

texts:
How can I make a graph of the noisiness of some audio durations of 100ms - 10s with 100% being white noise and 0% being a sine wave.

I have FFT at the moment, are there other ways to know the commensurate random nature and perhaps even the dissonance of sound?


texts: What is the basic idea behind Fourier transform?

texts:
What is the basic idea behind (discrete and continuous) Fourier transform (FT)? In short, what is the difference between discrete and continuous FT?

I have read multiple answers on the web related to this topic, but they are all long, provide long examples or cryptic equations. I would prefer not to read another long answer or an answer with the example of the circles or ellipses. I do not care about that. Ideally, I would like to read short answers which give the idea behind this concept. Your target audience should be high-school students (or, definitely, not math or physics college students), or, in general, people that do not have a very solid mathematical background and do not want to read equations. Hence, I hope it is clear that equations should be avoided as much as possible, unless they can be understood easily by high-school students.


texts: Reason of effect of continuity and discreteness on periodicity of signal?

texts:
We know that a sine signal in continuous time is a periodic signal with a period $2\pi$ whereas the same sine signal in discrete time is aperiodic.
My question is how changing just a type of signal (by sampling) can affect the signal's periodicity and why this happens whereas we have just omitted some real values between samples and how can we make use of this fact practically?


texts: Convolution Theorem: Hamming Window on a Time Series and Fourier Domain

texts:
If we have a set of time series data, y, consisting of 100 data points. One can apply a N (odd) Hamming window as a weighted moving average to decrease the noise. Say, if we choose 7 point Hamming window, H,  [0.0800    0.3100    0.7700    1.0000    0.7700    0.3100    0.0800], and perform a weighted moving average on y. In Matlab, I am doing conv(y, H./sum(H), 'same').

Since this is a convolution, I am wondering that by convolution theorem, one should be able to obtain identical results in Fourier domain, because convolution should become multiplication in the Fourier domain.

However, the FFT of y will have 100 complex numbers, and Hamming window will also have 7 complex numbers. What would be an equivalent operation of a weighted moving average in the Fourier domain?

Thanks.


texts: Kalman filter for random observation matrix $G_t$

texts:
I have a problem that is similar to the state space model in Kalman filter but the observation matrix $G_t$ of $$y_t=G_tx_t +w_t,$$ is random. The elements of $G_t$ are i.i.d. random variables with a given distribution. Is it possible to use the Kalman filter to this problem? If it has already been used please provide a reference.

Above $w_t$ is noise vector of Gaussian distribution with known covariance matrix, $x_t$ is the state of the system at time $t$ and $y_t$ is the observation at time $t$. The state tarnsision is given by $x_{t+1}=F_tx_t+B_tu_t+v_t$ where $F_t$ is a known matrix and $v_t$ is Gaussian noise with known covariance, $B_t$ is also a known matrix and $u_t$ is Gaussian distributed random input.

Thanks a lot.


texts: Approximating pixel location in scale space

texts:
I'm implementing a scale-space with o octaves and s scale levels in each octave.
Each octave is half the size the previous.
I have keypoints found in each octave, I need to approximate the real location of the keypoint to their location in the original image (before the reduction). I need it to be more accurate then just doubling the coordinates.
I know how to approximate sub-pixel in each octave, but I need "between" the octaves.

Thanks
matlabit


texts: What is the range for a demodulated AM signal?

texts: Lets say I perform amplitude modulation (modulation index = 1) of a signal that has a range of [0, 1]. If I demodulate by taking the absolute value of this signal and then low-pass filter it, what will be the range of the demodulated signal (assuming the low-pass filter doesn't attenuate any of the modulation signal's frequencies)?


texts:
I'd like to generate a simple QPSK signal in GNU radio, but so far I haven't been able to get it to perform the way I'd like and I'm not sure if this has to do with some misunderstanding I have about how QPSK works or if I just am not familiar enough with GNU Radio.

Essentially, my conception of QPSK modulation is this - 2 bits are fed into the modulator (4 possible symbols - 00, 01, 10, 11) and, based on which combination of those 2 bits was received, one of 4 possible waveforms is generated (corresponding to +/-.707 +/- .707j). Thus, when using a QPSK modulator in GNU Radio I would expect to see, in the constellation diagram, a square with 4 points at each corner, but instead I see diamonds with thousands of points, even in the absence of any introduced noise. When I print the output of the QPSK modulator to a file sink, I also see more than 4 distinct points.

Is there a way to take a random input stream of 0's and 1's and produce a modulated signal from it, like QPSK, that has only one of four values? It is important for the simulation that I'm running that there are only 4 possible modulation values before noise is added to the system. Also, I will be sending these signals between USRP's, so if there is some way to generate a QPSK signal and transmit it on a carrier frequency between USRP's that doesn't rely on the PSK Mod block (I'm not too familiar with how to implement a raised cosine filter) that would be appreciated.


texts:
How can i get lpc coefficient by htk?
I know that there is a binary code "HCopy" which is ganna help with that... but my question is how can i read a wav file in htk?
and i guess i need activeperl to run this? right?
i installed htk but do i need compilation for getting lpc coefficient?
I would appreciate if some one introduce me a pdf or step by step instructions for doing this...
Thanks a million :)


texts: how to eliminate interference in the same center frequency

texts: I'm trying to receive an ASK signal in the 433 MHz center frequency with a bandwidth of 100KHz. I'm using a CC1101 transceiver connected to a Raspberry Pi GPIO ports. However, there is an intermittent signal in the 433MHz center frequency coming from the sensors of my alarm system, so the RX doesn't get all the samples of my wanted signal and I'm not sure how I could filter the interference because both signals are "unknown" (they are not generated by me). I wonder how I could do this, at least in theory.


texts: Taking FFT of 10 concatenated time traces with a time difference between each trace?

texts:
My goal - I'm interested in the white noise spectral region at higher frequencies, especially the phase information.
In this data acquisition instrument, digitised time data transfer is much faster for processing when "segmented", but with a twist that I am not sure about.
Example: Sampling rate of 60 MHz and ~1M samples acquired each segment
10 time-trace segments each with ~1 MSa (nearly 21.8 ms) but a "dead-time" between each segment of 10 us + small jitter up to 10 ns.
i.e. 10 segments of ~1 MSa with 10 us delay between each segment, total of 10 MSamples (~ 170.6 ms)
In post-processing I concatenate these 10 segments as one time trace of 10 MSa. But what is the effect of the 10 us delay between each segment in the frequency domain? I expect to see something along the lines of discontinuities... but what I'm be measuring is just white noise at high frequencies.
Alternatively: Just acquire one time trace of 10 MSa and FFT that.
Still learning Fourier transforms as I'm doing this and I'd like some advice please.


texts: How do I evaluate the performance of a segmentation algorithm?

texts: I'm currently on a project of liver tumor segmentation. I segmented the liver using region growing, and I have to assess how accurate the result is.  I recently learnt that there are certain metrics for assessing the accuracy of the region growing algorithm segmentation, like Tanimoto Coefficient, Correlation, etc. But I don't know how to implement them in Matlab. Check out https://stackoverflow.com/questions/9553204/tanimoto-coefficient-using-matlab


texts: Question based on scaling the FFT for obtaining PSD

texts:
Following is the code given in MATLAB's site to estimate PSD using FFT.

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));

N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

plot(freq,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)'

I know I might lack a proper understanding of the fundamentals, but can anyone explain in order to extract the first half of xdft why does the index run from $1$ to $N/2 + 1$ and not $1$ to $N/2$? (I suppose the Nyquist frequency lies at $i = N/2$, am I right?)
Some sources mention that the square of the magnitude should be scaled by $\frac{1}{N}$, while here it is scaled by $\frac{1}{Fs*N}$. I am unable to figure out which of these two should be used.
I understand that the DC value should be left untouched when converting to single-sided spectrum. But why is the scaling by $2$ performed only till end-1 and not till end?


texts: The origin of the restricted isometry property (RIP)

texts:
I have been looking for the origin of the restricted isometry property (RIP). Many papers cite the origin of the RIP in the following paper

Emmanuel Candès, Terence Tao, Decoding by Linear Programming, February 2005.

but I cannot understand where the definition came from. Some others papers say that the RIP is very similar to the Johnson-Lindenstrauss lemma (JL), and in fact it is, but I am not sure that that it is the origin of the RIP. However, I haven't been able to find this origin, mathematically.

On the other hand, how can the RIP be proved with certain matrices? Let's say for example Gaussian matrices or a Fourier matrix?


texts:
I want to calculate the average of a power spectrum in the complete frequency range. I have done a stft and want to calculate power in dB and then the average. Please let me know from the code below if my understanding is correct and how to calculate the average of the power spectrum.

y, sr = librosa.load('Noise Recorder.wav',sr=None,duration=2)
n_fft = 8192
S = librosa.stft(y=y,n_fft=n_fft,win_length=n_fft,hop_length=n_fft,window='hann')     # Short FFT
F = librosa.fft_frequencies(sr=sr, n_fft=n_fft) # Frequency components
Abs = np.abs(S)  # Same as mag, phase = librosa.magphase(S)
Abs_avg = np.mean(Abs, axis=1)
Power_avg =librosa.power_to_db(Abs_mean**2) # dB
Scale_dB = 86.85    # Scaling factor

Power_Scaled_dB = Power_avg + Scale_dB  #Plot #Max:114dB at 1KHz
RMS = np.sqrt(np.sum(Power_Scaled_dB**2) / n_fft) # n_fft=8192 #Ans:24.4

###### Additional RMS tryouts #####
# Method 1 : Based on PSD formula
Px = 10**(np.max(114/10)) # Ans: 251188643150.9
RMS_Px = np.sqrt(np.sum(Px**2) / n_fft) # Ans: 2775268639.5
RMS_Px_dB = 10*np.log10(RMS_Px) # Ans: 94

# Method 2: Based on amplitude fomrula
A = np.sqrt(Px) # Ans: 501187.2
RMS_A = A / np.sqrt(n_fft) # nfft= 8192(fft length), Ans: 5537.388930343294
RMS_A_dB = 20*np.log10(RMS_A) # Ans: 74

With a 114dB 1Hz Sine wave input signal, after stft i obtain the psd in the plot
Now I want to derive the PSD from the spectrum, to further calculate the average/RMS of the spectrum. Given by the relation PSD and from there I want to calculate RMS given by the relation RMS/Leq. I am unsure how to derive the PSD from the Power_avg
 and


texts: Eigen Value and Eigen Vectors in Image Processing

texts: I have been reading about eigen values and eigen vectors but i haven't be able to find any decent explanation relating their application in image processing/computer Vision for example say face detection and eye detection. Can somebody please through some light on this?


texts:
The wikipedia entry on PSD has one confusing line:

  Summation or integration of the spectral components yields the total power (for a physical process) or variance (in a statistical process)

But would the total power correspond to the uncentralized second moment, so variance + mean2? In particular, there are stationary processes with a nonzero mean. Is the integral over the PSD still just the variance?


texts: Why can't we have IoT from 2005? but begin now,or future?

texts:
In 2005,ITU announce the "IoT",but why can't we have IoT at that time?

Is it because the speed of internet is not fast enough?wireless communication is not good at that time?Does anyone have some good reason which can explain that we just have the IoT concept in 2005,but we can not have the IoT in the real world at that time?


texts:
My understanding of what an OFDM modulator does, is that it takes a parallel sequence of QAM symbols of size M, pads it with zeroes corresponding to the zero subcarriers, then passes them through an IFFT block of same size M. IFFT is used because it can digitally create subcarriers.

My questions

1- Is the result output of the IFFT operation, an OFDM symbol of length M too?

2- Is EACH OFDM symbol going to be weighted sum of the ALL the QAM symbol?

3- How do things change if we have two antennas instead of one, do we have an OFDM modulator to each antenna? For example when using the Alamouti scheme? I have no clear vision how things change in this case, i.e what is the input at the IFFT now?

I am looking forward for answers thank you.


texts: Difference between recursive and system with memory

texts:
What's the difference between a recursive system and a system with memory?

$$y[n] = y[n - 1] + x[n]$$

Is a recursive and with memory, can a system have no memory and recursive?


texts:
I have a square-ware signal with 1V peak-to-peak (alternating between -500 and +500mV, so its absolute value is a constant 500mV) and I would like to add to this an AWGN expressed in dBc ("decibels relative to the carrier").

How I should generate such normal distribution noise?

I got a hint that - for e.g. -40dBc - the following code would work:

normal(0, (m.pow(10, (-40 / 20)) / 2), ..)

However this ends up generating a lower dBc value (around -31 dBc), meaning "stronger" noise than expected.

Reading around made me find the description of "half-normal distribution", which states that the expected absolute value scales not with sigma, but with sigma*sqrt(2/pi). I tried the following code - and it does seem to work - but I would like to ask for your kind input on whether my thinking was correct and if this code is the appropriate one:

normal(0, (m.pow(10, (-40 / 20)) / 2 / m.sqrt(2 / m.pi)), ..)

The way I calculated effective dBc was this:

Integrate carrier (in this case it is t*0.5 [V] <- let us call this a
Integrate the output of normal() <- let us call this b
Calculate 20 * log10(b / a)

Notes:

The code sniplets above are Python-like
The "/ 2" is there as my carrier is 1V peak-to-peak


texts: Are there any cases where the negative frequencies of an FFT are used?

texts: In performing FFTs and looking at spectra for real signals (say, taken from a sensor) we ignore the negative frequencies. I understand that the negative frequencies often reflect the positive frequencies because their imaginary components are complex conjugates. Is this only the case for real signals? Would we be able to notice a difference if imaginary components of a signal did not cancel out (as in hearing an audio signal)?


texts: sinus generation and it's fft in radians

texts:
I'm trying to generate sine signal and display it's fft in radians, byt I meet several issues that I can't understand. Please find my script below:

clear all; close all; clc;
fs = 10*pi;
ts = 1/fs;
t = 0:ts:2*pi;
sine = cos(10*t);
plot(t/pi,sine)
sinefft = fft(sine);
subplot(2,1,1)
plot(abs(sinefft))
subplot(2,1,2)
plot(abs(fftshift(sinefft)))

Results

1. Why on first plot fft shows that cosine has a frequency 11 radians?

2. How to shift fft to achieve frequency domain around 0?

I will be thankful for your explanations.

Maks


texts: Correct way to simulate Rayleigh fading for indoor environments

texts:
I am trying to simulate indoor communications with all components included. But something disturbs me when it comes to Rayleigh fading part.

I am adding my rayleigh fading to the transmitted signal with that code:

rayleigh_h = 1/sqrt(2)*[randn(1, length(s_l) / OSR) + 1j * randn(1, length(s_l) / OSR)]; %Constructing rayleigh channel complex coefficients
rayleigh_h = repelem(rayleigh_h, OSR); %Coefficient repeating to apply rayleigh fading to all samples
s_l = s_l .* rayleigh_h;

s_l: the complex modulated signal

OSR: Over sampling ratio

If I need to clarify the code a bit: Every bit is being sampled by a number of samples named "Over sampling ratio". Since I want to add rayleigh fading symbol-by-symbol instead of sample-by-sample, all samples that belong to a symbol are multiplied by the rayleigh coefficient corresponding each bit.

Sorry for if the code clarification is a bit confusing.

Here is my real question.

The method I am using above obviously changing the amplitudes of the signal with multiplying the with complex gaussian coefficients at instants. However, at indoor locations, it is fairly correct to assume that mobility(speed) of the devices should be limited with the walking speed. Therefore changes with the signal, and the fading coefficients should change slowly, not that fast. I am trying to figure out how to simulate it.

Any help is appreciated. Please also advise me if I have misunderstood the concept of the rayleigh fading.


texts: What is the maximal frequency resolution for Matlab's STFT implementation spectrogram()?

texts:
Matlab's spectrogram() function calculates the STFT of a signal. It describes its NFFT argument as follows:

  S = SPECTROGRAM(X,WINDOW,NOVERLAP,NFFT) specifies the number of
     frequency points used to calculate the discrete Fourier transforms.
     If NFFT is not specified, the default NFFT is used.

Am I correct in that NFFT is a trade-off only between frequency resolution and number of computations? For my offline work, there's no need to save cycles. Is there any maximum limit for NFFT, imposed e.g. by spectral leakage, or any other problem that I should know about, or can I set that argument to as high as possible?


texts: What is the effect of reducing NFFT less than the signal length

texts:
I am analyzing a signal in MATLAB by taking its fft. The singal has $200$ samples and I am currently taking nfft = 256. I understand that this would result in zero padding. Now, I want to try to reduce the nfft to 128, 64 etc. I wanted to ask what will happen to the fft.

Matlab says for Y = fft(X,n), If X is a vector and the length of X is greater than n, then X is truncated to length n. What would this actually result in. Can anyone explain, please?


texts:
Determine whether the following system is time-invariant or not:

$y(t)=x(t)\sin 10\pi t$

Solution:

Given: $y(t)=x(t)\sin 10\pi t$

$y(t)=T[x(t)]=x(t)\sin 10\pi t$

The output due to input delayed by $T$ sec is:

$y(t,T)=T[x(t-T)]=y(t)|_{x(t)=x(t-T)}=x(t-T)\sin 10\pi t$

The output delayed by T sec is:

$y(t-T)=y(t)|_{t=t-T}=x(t-T)\sin 10\pi (t-T)$

$y(t,T) \ne y(t-T)$

Conclusion: The delayed output is not equal to the output due to delayed input. Therefore, the system is time invariant.

  My query/doubt:

When the output delayed by T sec :

$y(t-T)=y(t)|_{t=t-T}=x(t-T)\sin 10\pi (t-T)
=x(t-T)\sin  (10\pi t-10\pi T)$

In case the value of $T$ is $\frac{2n}{10}$; $n\ge 0$, then the above equation becomes equal to $y(t,T)$

ie,

$y(t,T) = y(t-T)$ , When $T=\frac{2n}{10}$; $n\ge 0$

So can we say that the system is time-invariant when $T=\frac{2n}{10}$; $n\ge 0$ ?


texts: FFT-based coarse carrier recovery for QPSK

texts:
An FFT-based coarse frequency acquisition block is being tested in GNU Radio. The algorithm used can be found on the GOES satellite user manual and is summarized in the figure below.

The algorithm works by:

Carrier detection by absolute-law, square-law or fourth-power law. Square law detector involves squaring the IQ signal while fourth-power law uses the fourth-power of the IQ signal.
Finding the FFT
Averaging successive FFT sequences
Finding the peak by looking at 8 adjacent bins with the highest magnitude.

I have been able to implement the algorithm in GNU Radio, at least for BPSK. The algorithm is able to recover the carrier at all offset scenarios (offset < 10% symbol rate, offset ~ symbol rate, offset > symbol rate).

According to the GOES manual, the same algorithm should work for OQPSK, and I assume by extension, QPSK. However, the carrier derived by the algorithm for the case of QPSK is pretty bad, especially in low-SNR scenarios.

I have two questions

Are there any literature sources as to how the frequency detectors (square, fourth-power) work? I have basically implemented the algorithms but I do not fully understand their mathematical derivations.
Is there a more effective frequency detector for QPSK apart from square/fourth-power/absolute laws that could remedy the situation?

Regards,
Moses.

EDIT

Disclaimer: Using a fouth-law as proposed by @DanBoschen resulted in expected BER performance. For all intents and purposes, the answer is already accepted. But, I have observed some behavior of the nature of FFT relating to this question. I figured there wouldn't be any need to open a new question since we are already here.

In the flowgraph, the sample rate used is 8Mhz. The fourth-law-derived carrier can be tracked from -4Mhz to +4Mhz. There seems to be an "FFT ambiguity" when the offset frequency is above $\frac{1}{8}$-th the sample rate i.e anything above 1Mhz. Taking, for example, an offset of 1.5Mhz. The fourth-law detector will create a carrier at 6Mhz, which will be wrapped to -2Mhz in the second Nyquist zone. Dividing this by 4 results into an offset of -0.5Mhz, instead of 1.5Mhz, as shown in the figure below (derived carrier in green).

We can solve this problem by simply using a sample rate which is 8x the maximum offset for the fouth-law or 4x the maximum offset for the square-law. But, of course, that's not an ideal solution. I was wondering if there are better ways.


texts:
I have a function which looks like so:

I would like to know if there is a measure for a "weak type run length".
How can I know that a function is approximately constant in some specific neighborhood?
How can I differentiate this function from random noise?


texts: Fitting VAR Process with Generalised Gaussian Noise

texts:
Consider the $m$-dimensional VAR process

$${\bf x}_t = \sum_{l=1}^{P} A_l{\bf x}_{t-l} + {\bf e}_t$$

where the componenets of ${\bf e}_t$ are spatially and temporally independent and follow a generalised Gaussian distribution i.e. $e_{i,t} \sim GGD(\alpha_i, \beta_i)$.

Given my data $X$, I am trying to estimate the VAR coefficients $A_l$ and GGD parameters $\alpha_i, \beta_i$. I have tried to use Newton-Raphson optimisation with the vector

$$[\alpha_1, ..., \alpha_M, \beta_1, ..., \beta_M, A_1(1,1), ..., A_1(1,M), A_1(2,1), ..., A_P(M,M)]^T$$

but it is painfully slow and not very robust.

What would be the best way to tackle this model fitting problem?


texts: Recursively obtaining all-pole coefficients from complex cepstral coefficients?

texts:
This is a problem from the book Discrete-Time Speech Signal Processing by Thomas Quatieri, Chapter 6 Problem 13:

  Suppose $h[n]$ is all-pole (minimum-phase) of order $p$. Show that the all-pole (predictor) coefficients can be obtained recursively from the first $p$ complex cepstral coefficients of $h[n]$.

I was only able to solve for the opposite. i.e. complex cepstral coefficients from the pole coefficients.


texts: Requirement on processor and memory for image processing

texts:
I'm designing a camera system to use a 5MP CMOS image sensor with a micro controller board. In order to select proper components I'm looking for some numbers. It would be very appreciated if you give me advices.

The question is this. If I take a picture with resolution of 2592 x 1944, frame rate of 1 fps and the format of image data is 32 bit RGB, then how fast cpu and how large memory would you need? When I calculated the bitrate, it's like the followings

data size of one frame
number of pixel = 2592 x 1944 = 5,038,848
data size of one frame = 32 bit x 5,038,848 = 161,243,136 bit = 157,464 byte
(And I confirmed this size with raspberry pi camera)
data transfer
Now I need to transfer the frame data to my flash memory. But I would take just one cut with long exposure time. In other words, I don't need to take 2 or 3 pictures for short time. It would be enough to take just 1 image for 1 hour. Let's say I take 1 picture for 1 seconds and don't take any pictures anymore. Then,

bitrate = 161,243,136 bit/s = 157,464 byte/s = 15.4MB/s

But I can't go ahead with this number. Does it mean that I need a RAM larger than 15.4MB to save the frame for a moment? Or is it possible to use 2MB to transfer the data to hard disk or some storage, even it's slow? I couldn't find any information by googling.

In addition, I wonder how fast processor I need. For example, do I need 200MHz processor? or is it enough to have 15MHz processor like arduino? I have no idea. Please help me.


texts: Overlap/Chunk processing possible with Zero-phase?

texts:
I'm relatively new to DSP and I'm currently investigating FIR filter and IIR filters. From what I've found FIR filters can be implemented efficiently using the overlap-save method, but I was wondering if it might be possible to do the same using a IIR zero-phase filter?

For example say I have a fairly long signal if I divided the signal into overlapping chunks similarly to how overlap-save does it, and I zero-phase filter each of these chunks to prevent phase distortions then discard the overlapping samples as per the algorithm, is there any way I would be able to end up with a consistent end result? I've been looking around the web and in books but I haven't been able to find any sources that says it's possible except for this site:

http://www.dsprelated.com/freebooks/sasp/Overlap_Add_OLA_STFT_Processing.html

Where it says for Zero-phase to overlap the segment with two of it's neighbors. I've tried experimenting with the instructions and attempted to implemented and used an overlapping low-pass zero-phase filter on a signal as seen below, with the overlap length being 150 on each end of the segments. I've also padded the borders of each segment with 0s to prevent border distortions. From the looks of the result I got it doesn't seem like that it'll be smooth however when I combine them compared to if I filtered it all in one go.

If anyone could give me some pointers on what I might be doing wrong and where I need to look further it'd be great thanks!


texts: Why is the DFT used for spectrograms rather than the DCT?

texts:
Related to Could a DCT be used for an audio magnitude spectrum rather than DFT?, but more on the (audio) spectrogram side of things

I've been reading about the DCT (such as from here) and it seemingly provides quite a few advantages over the DFT such as having more frequency bins for the same number of samples and better power compaction, which sounds useful for signal analysis.

Yet all the programs I've used that have spectrograms seems to use the (magnitude of) DFT rather than DCT. Why is that? Is there something bad or "missing" from the DCT that is important for spectrograms?

(Note: I already understand the benefits of the DCT for compression, and this question is about audio (or other time-displacement) signals, so please avoid focusing too much on compression or images)


texts:
I am working on a spectrum analyzing simulator, which is meant to compute and plot the FFT of a sine wave carrier signal of $16 Hz$.

From my practicals with a spectrum analyzer, I know that the spike shows up on the $16Hz$ mark. But after computing my own FFT and passing the magnitude of the FFT to a plotter in Matlab, I get two spikes, one at the beginning and the other at the end.

My sample rate is $2048$ samples per second, and what I plotted was the respective $abs(FFT)$ against a frequency step of $16/2048$. Doing the fftshift only moves the spike to the $8Hz$ mark and not $16Hz$.

I've tried computing FFT in different programming languages and with different tools (java jtransform, Matlab, C language, Excel, two java written classes for FFT by different universities).

I dont think my FFT result are wrong, I think it has to do with the plotting. I do not understand why my plot does not look like it would in real life. From my description and information, can you see any problems that would cause my plot to be wrong?


texts:
I am trying to shift an image on a subpixel level, with the method described here:
http://www.mathworks.com/matlabcentral/fileexchange/18401-efficient-subpixel-image-registration-by-cross-correlation/content/html/efficient_subpixel_registration.html

The results are however disappointing. Since I am actually programming this algorithm for the rotation of images, where I resample inbetween pixels, this problem keeps on bugging me: http://spin.ecn.purdue.edu/fmri/PDFLibrary/TongR_MRM_1999_41_253_256.pdf

Can somebody tell me how these strange artifacts can be resolved? According to multiple scientific articles I found, it should be possible to resample on other grids with satisfactory results. My minimal code example (in python) is below, along with the resulting images, which are not satisfactory at all.

Does someone know what causes this effect, and how it can be resolved to reproduce the results others have?

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

N = 2048

img=mpimg.imread('koe.png') #http://www.clipartbest.com/cliparts/RiA/yXL/RiAyXLriL.png

f = np.zeros((N,N))
f[:img.shape[0],:img.shape[1]] = img[:,:,3]

def fft(x):
    return np.fft.fft2(x)

def ifft(x):
    return np.fft.ifft2(x)

exp = np.exp
pi = np.pi
j = 1j

def s(r):
    plt.figure()
    #plt.imshow(np.real(r), vmin=0.0, vmax=1.0, interpolation='none')
    plt.imshow(np.absolute(r), interpolation='none', vmin=0.0, vmax=1.0)
    plt.colorbar()

x = np.tile(np.arange(N),(N,1)).T
y = x.T
result = ifft( fft(f)* np.exp(1j*2*np.pi*(-0.5*x/N-0.5*y/N)) )
s(result)
s(f)
plt.show()


texts: How to ensure anti-aliasing in image resampling at lower resolution?

texts:
Suppose I have a 1D image with intensities $f_i$, $i=[0,n-1]$, where $n$ is the number of pixels. I want to downsample it to have $\frac{n}{2}$ pixels (let $n$ be an even number). For each new point I want to average the intensity over 2 neighboring pixels as

$$\tilde f_i = 0.5(f_{2i}+f_{2i+1}), \qquad i=[0,n/2-1]$$

Can somebody, please, answer the following questions

Will such transformation introduce aliasing (and if so, what is the correct way to resample with averaging over pixels?)
Is there some math procedure, using which I can answer questions like #1 myself?
How, in general, discrete signals are anti-aliased? From DSP books I know that for continuous signals you just apply low-pass filter to your input signal and then sample whatever you got after the low-pass transform. But here is the problem: for images from the very beginning I have only discrete signal and I cannot figure out how to use continuous coordinate low-pass filters for discrete signals. May be I should treat an image as a piece-wise continuous function, which is constant within each interval $[i,i+1)$ ?


texts: Matlab: White noise with flat/constant power spectrum

texts:
I am relatively new to signal processing techniques and Matlab and need a bunch of test data in the form of white noise, as defined on the Wikipedia page, with constant, flat power spectrum.

I tried the matlab built-ins wgn(m,n,p) and randn(m,n) to get vectors of white noise, and then used different combinations of the periodogram/pwelch/dspdata.psd functions to check what I believe to be the power spectrum. I cannot say I am 100% certain on the basis of documentation and other forum sourcing, so do correct me if I am mistaken.

The resulting plots I get do not look even close to a flat spectrum but instead inherit the "zick zacky" curve of the original random data. The question I inevitably had to ask myself: Do you even know what to expect and how to understand the results of the functions and how a flat power spectrum looks like? I don't, hence I am here, looking for clarification.

I am looking for answers which go beyond the Wikipedia pages etc., as they left me confused.


texts:
In  the  big  data  era,  in  order  to  control  the  cost,  complexity,  and  bandwidth  of  collecting  and  processing  high-dimensional data systems, it is critical to exploit models that encapsulate prior information regarding the signals of interest
The compresive sensing is a new signal processing technique which may use it instead of nyquist sampling theoreme, it's about sensing and compressing in the same time.

Compresed sensing has many applications indcluding MRI, RADARS,..this new technique rely on sparsity of the signal of interest.

In addition of sensing and compressing signal we need to reconstrut the signal when we would like to apply some process.., for this reason there are several algorithms that can be classified into 3 types:
*Optimization Methodes
*Greedy Methods
*Thresholding-Based Methods
Till here it's okay

my questions are :
What is(are) the most efficient recovery algorithms ?
What are the metrics used to define performance of a recovery algorithm ?


texts: What is the inverse Z transform of this:

texts:
$X(z) = \displaystyle \frac{1}{z}{\left(1-\frac{z^2}{4}\right)\left(1+\frac{1}{z}\right)\left(1-z\right)}$

Using partial fractions expansion i came up to this:

$\displaystyle \frac{1}{X(z)} = \frac{-\frac{z}{3}}{1-\frac{z}{2}} + \frac{\frac{z}{3}}{1+\frac{z}{2}} + \frac{\frac{2z}{3}}{1-z}$

I tried but i can't go on any further.


texts: Finding the fundamental frequency of a periodic signal

texts:
Suppose we have the signal $$x(t) = e^{j\omega_1 t} + e^{j\omega_2 t} + e^{j\omega_3 t},$$ where all the frequencies are rationally related (that is, the ratio of any pair of frequencies is a rational number).

How do I prove the fact that the fundamental frequency of the expression above is given by $\omega_0 = \gcd(\omega_1, \omega_2, \omega_3)$?

My attempt is as follows: assume that the signal has a fundamental period $T_0$ which is related to the fundamental frequency by $\omega_0 = 2\pi/T_0$. Now for the signal to be periodic with this period as the fundamental period the following equations must be satisfied $$\omega_1 T_0 = 2 \pi l$$ $$\omega_2 T_0 = 2\pi m$$ $$\omega_3 T_0 = 2\pi n,$$ where $l, m, n \in \mathbb{Z}$ and their values are such that $T_0$ is minimized.

Where do I go from here?
How do I show that the fundamental frequency $\omega_0 = 2\pi/T_0$ is the greatest common divisor of $\omega_1$, $\omega_2$, and $\omega_3$?


texts: Discrete-time convolution response on matlab using for loop

texts:
An LTI system's Unit impulse response is defined as:

> n = 1:100; h = 0.08*(0.9).^n;

if we give input signal as:

> x=[zeros(1,100) ones(1,100) zeros(1,100) ones(1,100) ones(1,100)];

to this system. What is going to be the output signal?
And how to find this using discrete-time convolution with a matlab code using for loops?


texts:
I'm curious about the algorithm and basic principles behind Optical Character Recognition.

Can anyone explain:

The most basic OCR algorithm that can be fairly easily implemented from scratch? Or the common denominator between different OCR algorithms?
Extensions/changes made to tailor the algorithm for application specific goals?

Also see my follow up question on OCR.


texts: Image super resolution algorithm in Matlab

texts:
I'm trying to implement a simple Image super resolution algorithm (DWT-Based Resolution Enhancement
) in the following paper

http://www.ripublication.com/aeee/52_pp%20%20%20405-412.pdf

I tried to implement the algorithm in figure 3 of this paper using Matlab.Code is given below.

img1 = imread('lena1.jpg'); %original High resolution image
[height, width, dim] = size(img1);

%%Downsampling the image by averaging
avgfilter = fspecial('average', [2 2]);
avgimg = filter2(avgfilter, img1);
img = avgimg(1:2:end,1:2:end); %Input low resolution image

[LL,LH,HL,HH] = dwt2(img,'haar'); %Decomposing

%Bicubic interpolation by factor 2 on each subbands
LL1 = imresize(LL,2,'bicubic');
LH1 = imresize(LH,2,'bicubic');
HL1 = imresize(HL,2,'bicubic');
HH1 = imresize(HH,2,'bicubic');

%% Calculating Difference image
for i=1:256
    for j=1:256
        img3(i,j,:) = img(i,j,:) - LL1(i,j,:);
    end
end

for i=1:256
    for j=1:256
        LH13(i,j,:) = img3(i,j,:) + LH1(i,j,:);
        HL13(i,j,:) = img3(i,j,:) + HL1(i,j,:);
        HH13(i,j,:) = img3(i,j,:) + HH1(i,j,:);
    end
 end

%bicubic interpolation(Here alpha = 2;Hence alpha/2 = 1)
 img31 = imresize(img3,1,'bicubic');
 LH131 = imresize(LH13,1,'bicubic');
 HL131 = imresize(HL13,1,'bicubic');
 HH131 = imresize(HH13,1,'bicubic');

img4 = idwt2(img31,LH131,HL131,HH131,'haar'); %IDWT
t = uint8(img4)
imshow(t);
imsave;

Input image that I used is given below

But I'm getting a completely unexpected output as follows.

Why this is happening.Please help.Thanks in advance.


texts:
I am following the book The Intuitive Guide to Fourier Analysis & Spectral Estimation with MATLAB. I am trying to selflearn the fourier analysis in matlab.
I got lost in one passage in the demonstration that states that the FT of the ACF function is the square of the DTFT of the signal. I have attached it here As you can see in the passage that I named 1. a delta_tau is missed in my opinion. Can you confirm that? More important in the point 2. in my opinion there should be a tau and not t. So I don't know how it is demonstrated this formula, because if in the passage 2. there is t the demonstration can be ended. I hope you can help me


texts: Why is the length of smoothed data from filtfilt always double of the input?

texts:
I am trying to filter the data I have using python with the code :

import scipy.signal as signal

#First, design the Buterworth filter
N = 3     # Filter order
Wn = 0.1  # Cutoff frequency
B, A = signal.butter(N, Wn, output='ba')

smooth_data = signal.filtfilt(B, A, np.ravel(rawdata))

plt.plot(rawdata, 'r-')
plt.plot(smooth_data, 'b-')
plt.show()

but the length of the smoothdata is always double the input ! so the time of the output signal will be doubled. which is not what I want.

is there any reason and solution for this ?


texts: Problem in constructing complicated block matrix

texts:
So, Here I want to solve this above equation and determine $(W|a_1 a_2 a_3)$ in Python.
However, most of the part I can possible do but I am totally clueless how could I construct such a Block Matrix ($L$) provided I have constructed all its constituting elements such as $P, K$ and transpose of $P$ also.
I was thinking of using reshape(), however could not guide my self to an ultimate solution.


texts: How to Apply Statistical Algorithms of Signal Processing to Regulate Variation of a Curve?

texts:
Below I am posting 2 graphs. I want to regulate the curvature of first graph using some statistical methods such as use of standard deviations, and modulate my graph to look like second one.

I am not asking here for any matlab command, but for some statistical method or a signal processing algorithm I can implement using MATLAB.

I just want transformation of shape as that of image 2, I do not want to achieve exact values of image 2.

Thanks & Regards, Nick.

.


texts: Any simple Digital All-Pass Filter Design Tools?

texts:
I would like to be able to implement an All-Pass filter with a frequency-dependent group delay.  I need a maximum group delay at low frequency of about 20 samples (for Fs = 44.1kHz), and it needs to fall to zero(-ish) at the Nyquist frequency.  Ideally, I would like to specify a corner frequency and a rate at which the delay falls off either side of it.  How the heck would I design such a thing?

I'm afraid I'm not sufficiently competent to follow the design approaches I can find on the web that plunge straight into poles, zeros, transfer functions and complex conjugates.  What I can do, though, is implement the darned thing efficiently in code once it is designed.

I am able to use Matlab to design FIR and IIR low-pass filters, because it has a built-in tool that spits them out.  But it doesn't have an All-Pass filter design tool (at least not one I can find).


texts: What do troughs represent in a sound wave?

texts: What's the notion of negative amplitude? I hope the title and the question don't mean different things.


texts: Continuous-time RNN and Shannon sampling theorem

texts:
The most-used discrete-time RNN equations used in Deep Learning these days are those of Elman:

I have seen two very different continuous-time version of these, with different justifications. The most "natural" to me are the leaky integrator onces (here sigmoid is replaced with tanh, but any activation function would do):

which is justified by saying that taylor expanding this equation with discretization step of 1 gives back the discrete time equations above (see this).

Another version is found on Wikipedia (notice how the weight is outside the activation function!)

and justified by saying

  Note that, by the Shannon sampling theorem, discrete time recurrent neural networks can be viewed as continuous-time recurrent neural networks where the differential equations have transformed into equivalent difference equations. This transformation can be thought of as occurring after the post-synaptic node activation functions $y_{i}(t)$ have been low-pass filtered but prior to sampling.

Could someone explain this?


texts:
You have an experimentally-measured HRIR and an audio signal. And you wish to implement filtering using scientific libraries in Python such as numpy, scipy and matplotlib(for plotting)

Question A) In order to implement an FIR causal filter using a time domain convolution which function from the SciPy library do you use?

Question B) and what are the steps/stages involved between getting the DT samples and filtering- list or flowchart this answer.
Will post codes and plots here after these instructions

shared details for now:

The audio signal/samples is coming from a Wav file (1-channel) with 20 sec (882000 samples -rate :44.1 kHz)
HRIR are MATLAB (.MAT) files from a public database.


texts: linear convolution toeplitz matrix vs circular convolution toeplitz matrix

texts:
I have an issue in understanding the difference between building the Toeplitz matrix when the convolution is linear and when it's circular. As I know that Toeplitz matrix $H$ can be built as following
H = toeplitz(h;zeros(N-L,1),h(1), zeros(1,N-1));

where h is the channel and L is the length of the channel, and N is the total length of symbol convoluted with the channel.
My question, if the convolution with the channel is circular, will the Toeplitz matrix still be built in the same way?


texts: The meaning of correlating "Time Aligned" pseudo-random or walsh code sequences

texts:
Correlation is from a  mathematical point of view a shift in time domain +multiplication + summation operation

I thought that if 2 codes are "time aligned" then correlation reduces to a dot product multiplication and summation, so no shift.

So suppose I have a sequence -1     1     1    -1, then if it is time aligned with itself , the result of correlation or dot product +summation  is 4 and nothing else. I have to say that I acquired this understanding from looking at a CDMA receiver where I always  see a multiplication when signal are being despread

I guess my interpretation or understanding of the term "Time Aligned" when it comes to correlation is not correct, right?


texts: Weighted Sum of Auto Correlation - Lower Bound

texts:
Given a vector $ v $ with elements $ {\left\{ {v}_{n} \right\}}_{n = - \infty}^{\infty} $ and denoting $ {v}_{n}^{\left( k \right)} = {v}_{n - k} $, namely, a shifted vector by $ k $ elements (Mind the vector is infinitely long).

How could one prove that there exist $ \alpha > 0 $ such that:

$$ \sum_{k = -\infty}^{\infty} {2}^{- \left| k \right| } \left \langle {v}^{\left( 0 \right)}, {v}^{\left( k \right)} \right \rangle \geq \alpha {\left\| v \right\|}^{2} $$

One could see it a weighted sum of the Auto Correlation Function of the vector $ v $.

Hence can be written as:

$$ \sum_{k = -\infty}^{\infty} {2}^{- \left| k \right| } \sum_{n = -\infty}^{\infty} {v}_{n} {v}_{n - k} \geq \alpha {\left\| v \right\|}^{2} $$

My understanding of the question is that there some bounderay on how slow the Auto Correlation decays (There is nothing special about $ {2}^{-\left| k \right|} $ there, it can be replaced by any $ \frac{1}{c} > 1 $).

Thank You.


texts: How to compute the power bands of an eeg signal using python?

texts:
so I have an eeg signal (edf format) that has 25 channels and 248832 entries, sampling frequency of 512Hz. I have to compute the frequency bands:
– Delta: 0.1-4Hz
– Theta: 4-8Hz
– Alpha: 8-12Hz
– Sigma: 12-16Hz
– Beta: 16-36Hz
– Gamma: >36Hz
and plot them accordingly. I am using Python for this with scipy, numpy, etc. and I should get to something like this:

Does anyone have any point-outs/ideeas/tutorials that could help me compute the bands and then get such a plot(probably a histogram)? Thanks!


texts: Rayleigh energy equation - applied to Fourier Transform of a sinusoid

texts:
I saw details indicating "total energy of a Fourier-transformable signal equals the total area under the curve of squared amplitude spectrum of this signal."

I have no issue with that. But I have found pretty much no discussions about the consequence of applying the Rayleigh energy equation to the Fourier transform of sinusoidal signal - which of course has infinite energy. There are certainly discussions indicating that the energy equation can only be applied to finite energy signals, and therefore is not applicable to the Fourier transform of a sinusoidal time domain signal.

So - what is the bottom line? Does this mean that if one attempts to apply the energy equation to the Fourier transform of a sinusoid, then the result will be a finite value or finite number - and it will be a meaningless number, right?

That is - when applied to finite energy signals, the energy equation returns meaningful results. But when applied to infinite energy signals like a sinusoid, then it returns a value, but it will be a meaningless result, right?


texts: how can I analize music as spotify/musixmath/shazam does?

texts:
I mean, does exists some state of the art algorithm able to extract music features, classify them and identify similar music or find lyrics?

I know that there are many topics or steps involved on my question, like digital signal processing, feature extraction, machine learning, pattern recognition, etc. I'm a computer science student and I would like to know if there is a methodology to identify music as it does musixmatc. For example, this program is able to find the lyrics of a song from a very short audio sample, even with some noise. As I know, [1] is the unique related paper that is available.

1 An Industrial-Strength Audio Search Algorithm


texts:
I have been working on image processing. I have an image and then add it a Gaussian noise with standard deviation=0.005. then using averaging filter to denoise it (I know it is not a good idea). then I got negative PSNR. what this means? I think PSNR is a positive measurement.
the code for calculating PSNR is as follow:

function psnr1=PSNR(I,J)

mse1=(double(I)-double(J)).^2;
MSE1m=(mse1(:,:,1,:)+mse1(:,:,2,:)+mse1(:,:,3,:))/3;
psnr1=10*log10(255^2/mean(mean(MSE1m)));

end

what is wrong?
the final results for my code is these curves:


texts: Wind Turbine Vibration Analysis: Which approach?

texts:
I'm analysing vibration data from 4 wind turbines (WTs). 8 different sensors are sampled at 25.6 kHz for 10 seconds once a day. I have data from around 400 days (intervals). The plot underneath is for one of the four WT in the time domain.

The frequency domain from FFT of the four turbines (Gearbox HSS) in the the range 0-12kHz:

Current plan and questions:

I want to look for faults in the gearbox, and I am suspecting that WT 4 (bottom right plot above) is most degraded.
I want to apply some kind of high frequency technique for looking for early stage bearing damages in the higher frequency range, since these are not captured by FFT.

What approaches should can I consider here?  (I have looked a bit at Hilbert transform, but I'm not exactly sure how that would help me.)
I am also considering performing some kind of enveloping/wavelet transform.
Should I be looking at lower frequencies as well?

This is not my main study field, and I want to try building a machine learning classifier with variables obtained from the signal processing.

Plots of the lower frequency (0-2000Hz) development over the 400 intervals for all 4 turbines:

Thanks!


texts: Discrete Fourier Transform of PMF and PDF

texts:
I got this explenation:

"It is well-known that an image’s histogram is essentially
the probability mass function (pmf) of the image (only
differing by a scalar). Multiplying each component of the
pmf by a correspondingly shifted unit impulse results in the probability density function (pdf). Obviously, in the context
of discrete Fourier transform (DFT), the unit impulses can be
ignored, implying that we can treat pmf and pdf
exchangeable. Thus, the pdf can be thought as the
normalized version of a histogram."

and I am trying to fully understand it.

It is clear to me that histogram of an image are probability mass functions. How can we treat is as probability density function since it is discrete signal, not continuous?
Is the only reason that DFT of both PMF and PDF would give us same result so it doesn't really matter in this scenario?

And how can we even get PDF from PMF?
I don't get this part: "multiplying each component of the
pmf by a correspondingly shifted unit impulse results in the probability density function"
and I cannot find any good example or explenation.


texts:
The particle filter is based on the state and observation model equations

$x_{t+1}=f_t(x_t, v_t)$

$y_t=h_t(x_t, u_t)$

The idea is to randomly generate some particles then propagate them through the equations, resample, and normalize to get an estimate of the state $x_t$.

My question is: we need to know the state and observation model equations so why do we need the particle filter? I hope I am missing something because it seems like if this knowledge is required then why bother with particle filter, can't we just do ML estimation?


texts: aliasing folding phenomena in frequency domain matlab

texts:
Hello i am trying to create the folding phenomena of undersampling in matlab,
When i undersample the sampling frequency is 135 less than the Nyquist frequency for 70 Hz signal
ou will see that it is shifted back by the amount of this new sampling frequency (105-70=35) Hz.
as shown in the matlab plot bellow.
I want to create the folding effect shown in diagram bellow (shifted copy).
Where did i go wrong.
This is my code which producedthe plot in the first post.
my code is shown as one line for some reason.

f1=10;
f2=30;
f3=70;
% twice the sampling rate
Fs=1.5*70; % sampling frequency is a bit above 2 times to get all the peaks.
Ts=1/Fs;
Tn=0:Ts:1;
fft_L=length(Tn);
y4_samples=10*sin(2*pi*f1*Tn)+10*sin(2*pi*f2*Tn)+10*sin(2*pi*f3*Tn);
%stem(Tn_new,y4_samples);
ff=fft(y4_samples);
ff1 = abs(ff/fft_L);
fft2 = ff1(1:floor(fft_L/2)+1);
fft2(2:end) = 2*fft2(2:end);
f = Fs*(0:fft_L/2)/fft_L;
plot(f, fft2)


texts: Why FIR uses often direct-form instead of transposed-form?

texts: Even though tranposed-form has less computational burden and memory requirements, why FIR uses mostly direct-form?


texts: Is there any reason that cross correlation would perform well or poorly on sparse binary arrays?

texts: I am using matlab xcorr to correlate simulated photon count data that has some Gaussian random noise set on top of it and it is working fine when the average value in the arrays is greater than one but when I split the arrays into enough bins such that most elements are zero (conserving the total photon count) and the average is much less than one photon per bin, suddenly the lags from xcorr get much larger (relative to the total length even). For example, at one resolution the lag may be 10 while if I increase resolution ten times, the lags become ~500, much more than ten times the previous value. When I do this without any noise added it works very well, and I am thinking maybe because when the values are mostly zeros and ones and the rare two, any differences such as the xth bin being a one rather than a zero is a large percentage difference in the value and this leads to a large difference in the sums (~integrals) that xcorr calculates. Does anyone have any more knowledge or thoughts on whether this is correct or not?